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Abstract 

A  boundary  element  method  is  applied  to  the  prediction  of  the  flow  around  propeller 
blades,  with  emphasis  at  the  tip  region.  The  presented  work  is  divided  into  three 
major  parts.  In  the  first  part,  a  new  panel  arrangement,  namely  the  FLow  Adapted 
Grid  (FLAG),  is  proposed.  This  grid  is  normal  to  the  blade  leading  edge  outline  and 
adjusted  to  the  force  free  wake  geometry  at  the  trailing  edge.  The  location  of  the  tip 
vortex  detachment  point  is  determined  by  an  iterative  method.  The  effectiveness  and 
robustness  of  the  flow  adapted  grid  are  demonstrated  through  numerical  validation 
tests  and  through  comparisons  with  existing  experiments.  The  flow  adapted  grid  is 
found  to  improve:  (a)  the  predicted  velocity  flow  field  and  the  pressure  distribution 
at  the  tip,  (b)  the  convergence  of  an  iterative  pressure  Kutta  condition,  and  (c)  the 
overall  numerical  performance  of  the  method  and  its  consistency  to  lifting  surface 
theory.  The  second  part  addresses  an  algorithm  for  predicting  the  three-dimensional 
vortex  sheet  roll-up.  A  higher  order  panel  method,  which  combines  a  hyperboloidal 
panel  geometry  with  a  biquadratic  dipole  distribution,  is  used  in  order  to  accurately 
model  the  highly  rolled-up  regions.  For  given  radial  circulation  distributions,  the 
predicted  wake  shapes  are  shown  to  be  convergent  and  consistent  to  those  predicted 
from  other  methods.  In  the  final  part  of  this  thesis,  the  flow  adapted  grid  and  the 
three-dimensional  wake  sheet  roll-up  algorithm  are  combined  in  order  to  estimate  the 
propeller  loading/trailing  wake  interaction.  Predicted  forces,  circulation  distributions 
,  and  tip  vortex  trajectories  are  shown  to  agree  well  to  those  measured  in  experiments. 
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Chapter  1 


Introduction 

In  most  marine  propeller  applications  cavitation  usually  occurs  first  in  the  strong 
vortical  region  in  the  vicinity  of  the  blade  tip.  Inception  at  the  tip  may  occur  either 
on  the  blade  or  in  the  vortex  core  (downstream  of  the  blade) .  Knowing  the  details  of 
the  flow  on  and/or  behind  the  tip  is  thus  crucial  in  determining  cavitation  inception. 
Tip  vortex  cavitation  has  been  one  of  the  loudest  underwater  noise  sources.  Due  to 
this,  tip  vortex  cavitation  has  been  of  great  interest  in  naval  propeller  design.  A  tip 
vortex  is  most  likely  to  cavitate  when  the  blade  tip  is  subject  to  off-design  conditions, 
due  to  spatial  non-uniformity  of  the  inflow  to  the  propeller,  or  due  to  flow  inclination. 
In  order  to  delay  tip  vortex  cavitation  inception,  designers  often  unload  the  circulation 
distribution  at  the  tip,  thus  by  sacrificing  on  propeller  efficiency.  Therefore,  in  the 
design  and  assessment  of  propulsors  for  naval  applications,  it  is  essential  to  accurately 
predict  and  thus  control  tip  vortex  cavitation  inception. 

1.1  Objectives 

The  objective  of  this  thesis  is  to  develop  a  robust  and  efficient  panel  method  applicable 
to  general  shape  lifting  surfaces,  including  propeller  blades,  with  full  wake  alignment 
and  wake  sheet  roll-up  in  three  dimensions.  The  approach  is  to  use  a  new  panel 
arrangement  on  the  propeller  blades  to  obtain  accurate  prediction  of  the  pressure 
distribution  at  the  blade  tip.  In  addition,  a  higher  order  panel  method  and  a  redis- 
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cretization  scheme  in  the  calculation  of  wake  sheet  roll-up  are  used  to  obtain  more 
accurate  and  smoothly  rolled-up  geometry  of  the  wake  sheet.  Emphasis  is  placed  on 
obtaining  accurate  prediction  of  the  tip  vortex  trajectory,  the  geometry  of  the  wake 
sheet,  and  the  pressure  distribution  at  the  tip.  It  is  expected  that  such  a  solution  will 
improve  the  prediction  of  tip  vortex  cavitation  inception. 


1.2  Previous  Research 

1.2.1  Tip  Vortex  Cavitation 

The  primary  methods  for  predicting  tip  vortex  cavitation  inception  are  variations 
of  the  method  by  McCormick  [51], [50].  He  proposes  a  semi-empirical  approach  in 
which  he  considers  that  the  minimum  pressure  depends  on  the  near-tip  loading  and 
identifies  the  important  role  of  the  foil  tip  boundary  layer.  He  then  postulates  a 
power  law  relation,  S  ~  Re~^,  between  the  boundary  layer  thickness  6  and  the  local 
Reynold’s  number  Re.  The  constant  factor  multiplying  this  scaling  law  as  well  as  the 
exponent  are  geometry  dependent  and  thus  the  method  can  provide  predictions  only 
for  geometries  which  are  close  to  those  that  have  been  tested  in  either  model  or  full 
scale.  More  recently,  from  LDV  measurements  in  the  vicinity  of  the  tip  vortex  for 
several  planar  wing  configurations,  at  different  tunnel  facilities  and  flow  conditions, 
Fruman  et  al.  [15]  have  shown  that  the  velocity  field  close  to  the  core  of  the  vortex 
depends  on  two  parameters;  (1)  the  strength  of  the  tip  vortex  along  its  trajectory 
and  (2)  the  radius  of  the  vortex  core.  Having  these  parameters,  the  value  of  the 
minimum  pressure  inside  the  vortex  core  and  its  location  along  the  vortex  can  be 
readily  determined.  The  corresponding  pressure  coefficient  should  then  be  equal  to 
the  negative  value  of  the  cavitation  number  at  inception.  Recent  attempts  by  Dupont 
and  Cerrutti  [9]  to  apply  Reynolds  Averaged  Navier-Stokes  solvers  to  the  tip  vortex 
flow  have  led  to  poor  predictions  of  the  minimum  pressure  in  the  vortex  core,  mainly 
due  to  gross  overprediction  of  the  size  of  the  vortex  viscous  core  .  On  the  other  hand, 
panel  methods  have  been  found  to  be  useful  in  determining  the  minimum  pressure 
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coefficient  at  the  blade  tip,  which  then  can  be  correlated  with  the  cavitation  number 
at  tip  vortex  cavitation  inception.  It  is  believed  that  a  robust  panel  method,  with 
full  wake  alignment  ,  will  not  only  improve  the  accuracy  of  the  predicted  pressure 
distributions  at  the  blade  tip,  but  also  will  provide  the  foundation  for  predicting  the 
tip  vortex  evolution  which  will  ultimately  lead  to  more  reliable  estimates  of  tip  vortex 
cavitation  inception. 

1.2.2  The  Panel  Method 

Panel  or  Boundary  Element  methods  (BEM)  have  been  applied  for  the  analysis  of 
propeller  flows.  This  method  is  well  advanced  and  has  been  extensively  described  by 
Hess  and  Valarezo  [20],  Lee  [42],  Kerwin  et  al.  [33]  and  Hoshino  [23].  These  methods 
employed  a  potential  or  a  velocity  based  formulation.  The  investigation  of  different 
panel  methods  by  Lee  [42]  showed  that  the  perturbation  potential  based  panel  method 
is  the  best  for  propeller  applications. 

A  perturbation  potential  based  Boundary  Element  Method,  including  the  presence 
of  the  hub  and  duct,  was  developed  at  MIT  by  Kerwin  et  al.  [33]  and  Lee  [42]. 
This  method  was  a  low-order  BEM  based  on  Green’s  formula  with  respect  to  the 
perturbation  potential.  The  method  discretized  the  propeller  surface  and  wake  with 
planar  quadrilateral  panels  and  constant  strength  sources  and  dipoles  were  distributed 
on  the  panels.  The  panel  method  has  been  applied  successfully  for  the  hydrodynamic 
analysis  of  marine  propellers  (Lee  [42]).  In  most  application,  the  BEM  has  been  found 
to  predict  spanwise  circulation  distributions  which  are  consistent  to  those  predicted 
from  Vortex  Lattice  Method (VLM)  (Hsin  [24]).  This  means  that  the  circulation 
distributions  predicted  by  BEM  for  lifting  surfaces  with  thickness,  smoothly  (often 
linearly  with  thickness)  extrapolate  to  the  circulation  distribution  predicted  by  VLM 
for  the  same  lifting  surface  with  zero  thickness.  When  this  method  is  applied  to  a 
wide  circular  tip  propeller,  inaccurate  solutions  are  obtained  and  an  iterative  pressure 
Kutta  condition  (Kerwin  et  al.  [33])  has  been  found  to  diverge.  Most  importantly, 
the  prediction  of  the  pressure  distribution  in  the  tip  region  has  been  found  to  be 
inaccurate.  The  correct  prediction  of  the  pressure  distribution  is  important  not  only 


for  the  thrust  and  torque  calculations  by  integration  of  the  pressure  distribution  on 
the  blade  surface,  but  also  for  the  cavitation  inception  estimation.  A  grid  oriented 
along  constant  radii,  namely  the  conventional  grid,  has  previously  been  employed 
by  the  BEM.  This  grid  arrangement  results  in  high  aspect  ratio  of  the  panels  at 
the  propeller  tip.  A  blade  orthogonal  grid  (BOG),  which  is  orthogonal  both  at  the 
leading  edge  and  the  trailing  edge,  was  developed  by  Hsin  et  al.  [25]  in  order  to 
improve  the  resolution  and  to  reduce  panel  distortion  and  aspect  ratio  near  the  tip. 
The  blade  orthogonal  grid  was  found  to  improve  the  convergence  of  the  computed 
pressure  distribution  at  the  tip,  and  consequently  also  the  convergence  of  the  iterative 
pressure  Kutta  condition  (Hsin  et  al.  [25]).  Nevertheless,  when  the  blade  orthogonal 
grid  was  applied  to  lifting  surfaces  and  propellers  with  wide  circular  tips,  it  was  found 
that  the  circulation  distribution  behaves  non-physically  at  the  tip.  To  avoid  this  non¬ 
physical  behavior  of  the  method  at  the  tip,  the  FLow  Adapted  Grid  was  introduced 
by  Kinnas  et  al.  [36],  which  is  orthogonal  at  the  leading  edge  and  aligned  with  the 
resulting  mean  flow  at  the  trailing  edge.  Pyo  and  Kinnas  [64]  showed  that  this  grid 
improved  the  convergence  of  IPK  condition  and  the  behavior  at  the  tip  of  propeller 
blades. 

The  Kutta  condition  may  be  enforced  numerically  via  the  Morino  condition  [59]. 
This  condition  requires  that  the  dipole  strength  at  the  trailing  edge  in  the  wake  to 
be  equal  to  the  difference  of  the  potentials  at  the  trailing  edge  panels  on  the  blade. 
Lee  [42]  incorporated  a  correction  to  the  Morino  condition  in  two  dimensions  by 
applying  the  Kutta  condition  at  the  exact  location  of  the  trailing  edge.  Application 
of  the  Morino  condition  in  three  dimensions  does  not  always  guarantee  equality  of 
the  resulting  pressures  at  the  trailing  edge.  An  iterative  pressure  Kutta  condition 
based  on  Newton-Rapson  method  was  introduced  by  Kerwin  et.  al  [33]  and  extended 
to  unsteady  propeller  flows  by  Hsin  [24],  and  Kinnas  and  Hsin  [35]. 

The  geometry  of  the  trailing  wake  behind  propeller  blades  was  first  calculated  by 
Cummings  [8],  who  determined  the  wake  geometry  by  using  a  two-dimensional  time 
domain  approach  in  which  the  body  effect  was  not  included.  In  the  BEM,  the  geom¬ 
etry  of  the  trailing  wake  sheet  is  calculated  via  the  lifting  surface  method  of  Greeley 
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and  Kerwin  [16]  in  an  indirect  way,  where  the  geometry  was  determined  from  the 
requirement  that  the  trailing  vorticity  at  the  trailing  edge  and  far  downstream  should 
be  aligned  with  the  local  flow.  Since  the  local  flow  depended  on  the  trailing  wake 
geometry  in  a  nonlinear  way,  an  iterative  procedure  was  employed.  In  the  method 
of  Greeley  and  Kerwin  [16]  the  effects  of  coupling  between  the  thickness  and  load¬ 
ing  and  the  wake  sheet  roll-up  were  ignored.  In  addition  the  radial  contraction  and 
the  ultimate  radius  of  the  trailing  wake  were  assumed  to  be  given  from  experiments. 
In  summary,  this  method  only  included  axial  and  circumferential  velocities  at  the 
trailing  edge  and  the  ultimate  wake,  and  thus  suppressed  the  effect  of  wake  sheet 
roll-up. 

1.2.3  Vortex  Sheet  Roll-Up 

In  the  past  there  have  been  a  large  number  of  attempts  to  model  the  vortex  sheet 
motion  by  replacing  the  continuous  vortex  sheet  with  a  finite  number  of  discrete  vor¬ 
tices  or  alternatively,  by  replacing  the  dipole  sheet  with  segments  carrying  a  piecewise 
constant  dipole  distribution. 

Rosenhead  [66]  was  the  first  to  attempt  this  approach  with  an  analysis  of  the 
nonlinear  Kelvin-Helmholtz  instability  in  a  two-dimensional  vortex  sheet  of  constant 
strength.  Westwater  [74]  first  applied  the  discrete  vortex  method  to  the  problem  of 
vortex  sheet  roll-up  behind  an  elliptically  loaded  wing.  Attempts  to  improve  West- 
water’s  results  by  increasing  the  number  of  vortices  representing  the  vortex  sheet 
have  not  been  successful.  It  appeared  that  this  approach  inevitably  led  to  chaotic 
motion  in  the  region  of  the  tip  vortex,  which  resulted  in  loss  of  the  identity  of  the 
vortex  sheet.  Different  approaches  have  been  attempted  to  desingularize  the  solution. 
Chorin  and  Bernard  [7]  and  Kuwahara  and  Takami  [40]  introduced  a  finite  core  model 
for  the  vortices,  in  which  the  velocity  remains  finite.  Moore  [56]  used  a  process  of 
amalgamation  in  which  the  vortices  are  combined  when  they  approach  each  other  too 
closely  or  when  they  have  to  represent  high  curvature  regions  in  a  spiraling  sheet. 
Maskew[46]  employed  the  sub-vortex  technique.  Despite  all  these  improvements,  the 
discrete  vortex  representation  resulted  inevitably  in  numerical  instabilities  when  the 
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number  of  vortices  was  increased  or  the  core  radius  was  decreased.  In  the  case  of 
wing  tip  vortices,  the  method  of  rediscretization  by  Fink  and  Soh  [14]  has  been  more 
successful  in  obtaining  smooth  vortex  sheet  behavior  over  long  periods  than  have 
been  reported  previously.  However,  recent  investigations  by  Baker  [4]  on  the  stability 
of  the  rediscretization  method  for  the  case  of  double-branched  spiraling  vortex  sheets 
demonstrated  that  this  method  eventually  ended  in  chaos  as  well.  As  an  alterna¬ 
tive  approach,  Baker  [3]  suggested  cloud-in-cell  method.  In  this  method,  the  velocity 
field  due  to  the  discrete  vortices  was  computed  by  solving  Poisson’s  equation  for  the 
stream  function  due  to  a  grid-dependent  region  of  distributed  vorticity  in  the  two  di¬ 
mensional  plane.  The  main  disadvantages  of  this  method  are  that  the  method  and  its 
fine-scale  behavior  are  sensitive  to  the  size  of  mesh  ,  the  surface  boundary  conditions, 
the  number  of  vortices  and  the  time-step,  as  noted  by  Murman  and  Stremel  [60]. 
From  those  attempts,  it  appeared  that  the  discrete  vortex  method  was  not  adequate 
to  compute  smooth  vortex  sheet  roll-up  reliably. 

On  the  other  hand,  the  use  of  dipole  distributions  for  two  and  three  dimensional 
attached  flows  is  well  advanced  and  has  been  extensively  described  in  the  past  by  Hess 
[19]  and  Johnson  [28].  In  its  most  general  form,  the  body  is  represented  by  sources  and 
dipoles  and  the  kinematic  boundary  condition  is  applied  on  the  body.  Across  the  wake 
sheet  the  normal  velocity  is  continuous  but  the  tangential  velocity  is  discontinuous 
by  a  “jump”  corresponding  to  the  strength  of  the  vorticity  sheet.  Thus,  the  wake 
sheet  can  also  be  expressed  via  a  distribution  of  dipoles.  A  complete  description  of 
this  method  was  presented  by  Maskew  [47].  Hoeijmakers  et  al.  [22]  developed  high 
order  panel  methods  based  on  the  slender  body  approximation.  In  general,  the  panel 
methods  produce  smoother  vortex  sheets  than  those  of  the  discrete  vortex  models  of 
Maskew  and  Rao  [48]  and  Fink  and  Soh  [13].  A  number  of  fully  three  dimensional  flow 
models  have  been  developed  for  the  prediction  of  wake  roll-up  in  order  to  overcome 
the  limitations  of  the  slender  body  approximation.  The  most  well  known  methods 
among  the  high  order  panel  methods  are  the  Boeing’s  LEV-Model  by  Johnson  et  al. 
[29]  and  the  VORSEP-Model  by  Hoeijmakers  [22].  The  comparison  with  experiments 
has  been  good  to  encouraging  as  far  as  the  overall  accuracy  of  the  position  of  the 
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vortices  is  concerned.  The  local  accuracy  depended  on  the  shape  of  the  wing,  angle  of 
attack,  number  of  panels  and  the  type  of  smoothing  applied.  The  calculated  pressure 
distributions  on  the  body  compared  reasonably  well  with  those  measured  by  Hummel 
and  Redeker[26].  The  only  disadvantage  of  these  methods  is  that  an  initial  geometry 
of  the  vortex  sheet  had  to  be  assumed  or  deduced  from  analytical  solutions. 

Recently,  Reynolds-averaged  Navier-Stokes  (RANS)  methods  have  been  developed 
with  an  assumption  that  the  flow  at  an  infinitesimally  small  distance  upstream  is  con¬ 
ically  similar,  which  might  be  a  valid  approximation  for  viscous  flows  that  vary  slowly 
in  the  streamwise  direction.  Lee  [43]  applied  this  method  to  predict  the  roll-up  of  the 
leading  edge  vortex  on  delta  wings.  Ega  et  al.  [11]  solved  the  parabolized  Navier- 
Stokes  equations  on  elliptic  wings.  Stern  et  al.  [71]  applied  their  RANS  solver  on 
propeller  blades.  These  methods  have  demonstrated  capability  to  simulate  qualita¬ 
tively  the  physics  of  viscous  dominated  processes  such  as  flows  inside  the  core  and 
vortex  breakdown.  However,  the  current  status  of  turbulence  and  transition  mod¬ 
eling  is  not  yet  adequate  so  that  quantitative  predictions  (especially  for  the  pressure 
distribution  inside  the  vortex  core)  based  on  RANS  are  feasible. 

There  are  numerous  other  applications  which  do  not  fit  conveniently  into  the 
above  categories.  McCune  and  Tavares  [53]  solved  the  two-dimensional  unsteady  and 
large  amplitude  delta  wing  motion  problem  with  leading  edge  separation.  McCracken 
and  Peskin  [52]  combined  a  finite-difference  method  with  the  vortex  blob  algorithm. 
McAlister  and  Carr  [49]  solved  unsteady  vortical  flows  with  dynamic  stall  around  an 
oscillating  wing. 

1.2.4  The  Present  Method 

In  the  present  method,  a  potential  based  boundary  element  method  is  applied  for 
the  analysis  of  propeller  flows.  The  flow  adapted  grid  (FLAG)  is  developed  in  order 
to  solve  the  convergence  problem  for  the  typical  grid  arrangements  and  to  improve 
the  tip  flow  behavior.  The  geometry  of  the  trailing  wake  is  decided  directly  from 
the  panel  method  without  any  assumption  and  the  vortex  sheet  roll-up  is  included  in 
the  FLAG.  In  order  to  model  the  wake  sheet  roll-up  in  three  dimensions,  bi-quadratic 
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strength  dipole  distributions  and  hyperboloidal  panel  geometry  are  used.  Throughout 
the  numerical  calculation,  rediscretization  is  applied  as  a  smoothing  scheme. 

The  consistency  and  convergence  of  the  results  from  the  numerical  method  are 
validated  first,  and  then  the  method  is  applied  to  several  planar  wing  and  propeller 
blade  geometries  and  the  numerical  results  are  compared  with  existing  experimental 
data. 


Chapter  2 

Boundary  Element  Method 


2.1  Formulation 

The  fundamentals  of  the  BEM  are  described  by  Kerwin  and  Kinnas  [33],  Lee  [42], 
and  Hsin  et  al.  [25]  and  only  a  brief  description  will  be  given  in  this  section.  The 
method  is  based  on  the  classical  Green’s  third  identity  (applied  on  the  body  surface 
Sb)-- 


where  the  Green’s  function  G  is  the  unit  strength  source  in  three  dimensions;  (j)  is  the 
perturbation  potential;  Sw  is  the  trailing  wake  surface  as  shown  in  Figure  2-1. 


The  BEM  implementation  involves: 

•  Constant  strength  dipole  and  source  panels  on  the  blade  and  in  the  wake. 

•  Hyperboloidal  panel  geometry  (critical  for  highly  twisted  body  geometries). 

•  An  Iterative  Pressure  Kutta  (IPK)  condition  which  determines  the  appropriate 
strength  A(?!>  in  the  wake  in  order  for  the  pressure  jump  across  the  trailing  edge 
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to  be  equal  to  zero  at  all  spanwise  locations. 

Two  panel  arrangements  have  been  used  in  the  past:  the  conventional  grid  by 
Kerwin  and  Kinnas  [33]  and  Lee  [42]  and  the  blade  orthogonal  grid  (BOG)  by  Hsin 
et  al.  [25]. 


Figure  2-1:  Global  coordinate  system  fixed  on  the  propeller  blade. 

The  former  of  the  two  grids  and  the  proposed  flow  adapted  grid  (FLAG)  will  be 
described  in  the  next  sections. 


2.2  The  conventional  grid 

The  conventional  grid  has  been  used  traditionally  for  vortex-lattice  applications  on 
3-D  wings  and  propeller  blades  by  Lan  [41]  and  Greeley  and  Kerwin  [16].  It  has 
also  been  called  the  “constant  radii”  grid.  The  panel  edges  are  located  along  the 
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intersections  of  the  blade  with  cylinders  concentric  with  the  axis  of  propeller  rotation. 
The  geometry  of  the  trailing  wake  is  found  by  following  the  Greeley  and  Kerwin  [16] 
procedure  : 

•  The  tip  vortex  trajectory  starts  on  the  blade  at  the  highest  radial  position. 

•  The  wake  gridlines  are  aligned  with  the  axial  and  tangential  velocities  in  the 
wake. 


Figure  2-2:  The  conventional  grid  on  a  propeller  blade  and  its  trailing  wake. 

•  The  radial  contraction  of  the  wake  at  the  tip  is  an  input  parameter,  given  from 
experimental  information  (usually  equal  to  30  degrees). 

•  The  radius  of  the  ultimate  wake  geometry  is  also  an  input  parameter,  given 
from  experimental  information  (usually  equal  to  0.83i?). 
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The  corresponding  panel  arrangement  for  one  of  the  three  propeller  blades  (with 
40  chordwise  and  20  spanwise  panels)  is  shown  in  Figure  2-2.  The  spanwise  circulation 
distribution  predicted  by  applying  the  BEM  on  the  described  grid  for  propeller  N4119 
[27]  is  shown  in  Figure  2-3.  In  this  particular  case  it  took  14  iterations  for  the  IPK 
condition  to  converge  {\/S.Cp\te  <  10“®). 


Figure  2-3:  Circulation  distribution  on  the  propeller  N4119;  J  =  0.833.  Predicted  by 
applying  the  BEM  on  the  conventional  grid;  before  and  after  applying  the  Iterative 
Pressure  Kutta  condition. 

The  circulation  distribution  “before”  the  IPK  corresponds  to  the  Morino  [59] 
Kutta  condition  in  which  the  dipole  strength  in  the  wake  is  taken  equal  to  the  dif¬ 
ference  of  the  potentials  at  the  panels  at  the  two  sides  of  the  trailing  edge.  This 
condition  has  been  found  to  produce  pressure  distributions  which  do  not  match  at 
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the  trailing  edge,  especially  in  the  vicinity  of  the  tip.  The  circulation  distribution  “af¬ 
ter”  the  IPK  corresponds  to  the  modified  wake  dipole  strength  which  ensures  pressure 
equality  at  the  trailing  edge.  This  panel  arrangement  has  been  successful  for  conven¬ 
tional  geometries.  However,  when  this  grid  is  applied  to  extreme  geometries  such  as  a 
highly  skewed  propeller  and  a  propeller  with  large  tip  chord,  inaccurate  and  divergent 
solutions  are  obtained. 


Figure  2-4:  Circulation  distribution  on  a  circular  wing  planform  hydrofoil;  a  =  5.73°, 
Modified  NACA66  thickness  distribution  with  [r/cl^ax  =  0.2.  Predicted  by  applying 
the  BEM  on  the  conventional  grid  ;  before  and  after  applying  the  IPK  condition. 


This  is  because  the  conventional  grid  arrangement  results  in  high  panel  aspect 
ratios  and  highly  skewed  and  twisted  panels  at  the  propeller  tip.  For  a  circular  wing, 
which  has  extremely  large  tip  chord,  the  circulation  distribution  is  show  in  Figure 


2-4.  Notice  the  large  difference  over  the  span  between  the  circulation  distributions 
before  and  after  IPK  condition.  Especially  near  the  tip  the  peculiar  behavior  of  the 
circulation  distribution  after  IPK  condition.  In  addition,  for  a  propeller  with  large  tip 
chord  ,  the  pressure  distributions  at  three  spanwise  locations  and  the  mean  velocity 
vectors  at  the  first  control  points  in  the  wake  along  the  trailing  edge  are  shown  in 
Figures  2-5  and  2-6,  respectively. 


Figure  2-5;  Pressure  coefficients  predicted  from  BEM  applied  on  the  conventional 
grid  (after  IPK  condition);  propeller  N4119,  J=0.833. 

The  velocity  vectors  are  computed  from  the  superposition  of  the  inflow  velocity, 
U in,  and  the  velocities  induced  by  all  dipoles  and  sources  on  the  blades  and  all  dipoles 
in  the  wakes  of  the  blades.  Notice  the  singular  behavior  of  the  pressures  and  the  mean 
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velocity  vectors  in  the  vicinity  of  the  tip. 

2.3  The  blade  orthogonal  grid 

To  improve  the  resolution  and  to  reduce  panel  distortion  near  the  tip,  the  blade 
orthogonal  grid  was  introduced  by  Hsin  et  al.  [25]. 


Figure  2-6:  Mean  velocity  vectors  along  the  trailing  edge  of  the  propeller  N4119  ; 
J  =  0.833.  Predicted  by  applying  the  BEM  on  the  conventional  grid. 

In  this  grid,  the  expanded  blade  outline  is  assumed  to  be  given  as  a  cubic  B-spline 
curve.  Then,  the  arclength  along  this  outline  is  expressed  as  a  function  of  parametric 
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variable,  which  is  zero  at  the  intersection  of  the  leading  edge  with  the  hub  and  reaches 
to  the  maximum  value  at  the  intersection  of  the  trailing  edge  with  the  hub.  With  the 
tip  point  defined,  the  outline  of  the  blade  is  divided.  Finally,  the  divided  points  are 
connected  by  cubic  B-spline  curves  with  four  vertices,  by  putting  the  second  and  third 
interior  vertices  along  the  normals  to  the  edge.  The  orthogonality  of  the  curves  can 
be  obtained.  Figure  2-7  shows  a  simple  example  of  this  grid  for  a  circular  planform 
wing.  One  can  see  that  the  grid  lines  are  orthogonal  to  the  outline  of  the  wing  and 
the  panel  aspect  ratios  remain  small  near  the  tip. 


Figure  2-7:  The  blade  orthogonal  grid  on  a  circular  planform  hydrofoil  and  its  trailing 
wake. 

When  the  BEM  was  applied  on  this  grid,  it  is  found  that  the  surface  pressures  at 
the  tips  of  non-lifting  bodies  were  computed  more  accurately  than  when  the  BEM  was 
applied  bn  the  conventional  grid.  This  is  the  consequence  of  concentrating  more  panels 
at  the  tip  as  well  as  of  producing  much  less  distorted  panels  (of  which  the  sides  are  of 
comparable  size  and  almost  orthogonal  to  each  other)  than  the  conventional  grid.  The 
blade  orthogonal  grid  is  also  found  to  improve  the  convergence  of  the  IPK  condition 
in  the  case  of  lifting  hydrofoils  or  propeller  blades.  This  is  a  direct  consequence 
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of  the  fact  that  the  trailing  edge  pressures  (which  drive  the  IPK  condition)  at  the 
tips  were  computed  more  accurately  now  than  in  the  case  of  the  conventional  grid. 
The  circulation  distributions  for  the  circular  planform  wing  are  shown  in  Figure  2-8. 
Notice  that  the  difference  between  the  circulation  distributions  before  and  after  the 
IPK  condition  is  now  larger  than  that  for  the  conventional  grid.  An  explanation  for 
this  will  be  given  in  Appendix  A. 


Figure  2-8:  Circulation  distribution  on  a  circular  planform  hydrofoil;  [r/c]max  =  0.2, 
a  =  5.73°.  Predicted  by  applying  the  BEM  on  the  blade  orthogonal  grid;  before  and 
after  applying  the  IPK  condition. 

Also  notice  that  the  circulation  distribution  after  the  IPK  condition  is  very  similar 
(also  “peculiar”)  to  that  in  the  case  of  the  conventional  grid.  The  results  shown  in 
this  and  the  previous  section  indicate  that  there  must  be  something  fundamentally 
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wrong  either  with  the  implementation  of  the  IPK  condition  and/ or  with  the  utilized 
grids. 

2,4  The  flow  adapted  grid  (FLAG) 

As  explained  in  Appendix  A,  the  difference  in  circulations  between  before  and  after 
IPK  condition  will  be  decreased  if  the  grid  on  the  blade  is  aligned  with  the  mean 
velocity  vector  at  the  trailing  edge.  This  means  that  if  the  grids  on  the  blade  and 
in  the  wake  are  aligned  with  the  mean  velocity  at  the  trailing  edge,  the  number 
of  iteration  for  the  IPK  condition  will  be  decreased  and  the  solution  will  converge 
after  fewer  iterations.  From  this  examination,  the  new  grid  arrangement,  namely 
the  “FLOW  Adapted  Grid  (FLAG)”,  is  developed.  The  flow  adapted  grid  will  be 
introduced  and  applied  on  three  dimensional  wings  and  propellers  in  this  section. 

2.4.1  Construction  of  the  grid 

The  flow  adapted  grid  is  constructed  by  satisfying  the  following  characteristics: 

®  The  grid  on  the  blade  is  adapted  to  the  resulting  flow  in  the  wake.  This  de¬ 
creases  the  number  of  iterations  for  the  IPK  condition  and  the  difference  in  the 
circulation  distributions  between  before  and  after  IPK  condition  as  explained 
in  Appendix  A. 

®  The  grid  on  the  blade  is  smoothly  connected  to  that  in  the  wake  at  the  trailing 
edge. 

@  The  grid  on  the  blade  is  orthogonal  at  the  leading  edge.  This  results  in  im¬ 
provement  in  resolution  near  the  tip  and  more  accurate  pressure  predictions  at 
the  leading  edge. 

@  The  grid  on  the  blade  includes  the  effect  of  the  resulting  flow  on  the  location  of 
the  tip  vortex  detachment  point  (also  called  the  “computational  tip”). 
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The  procedure  for  implementing  the  flow  adapted  grid  on  propeller  blades  consists 
of  the  following  steps  (also  shown  in  the  flow  diagram  of  Figure  2-9). 


Figure  2-9:  Flow  diagram  for  construction  of  FLAG. 


STEP  1 

•  Solve  the  boundary  value  problem  by  applying  the  panel  method  with  the  con¬ 
ventional  grid  and  straight  wake  for  the  wings  and  purely  helicoidal  wake  for 


the  propellers  (no  contraction).  Determine  the  dipole  and/or  source  strengths 
on  the  blade,  hub  and  wake  panels. 


®  Compute  the  total  mean  velocity  vectors,  at  the  first  control  points  in 

the  wake  at  each  radius,  via  the  following  superposition,  also  shown  in  Figure 


Figure  2-10:  The  contraction  angle  of  the  blade  gridlines  along  the  trailing  edge. 

V^AKE  =  Vt  +  Va  +  Vr-¥  Uir,  (2.2) 

where  V t,  Va,  Vr  are  the  tangential(circumferential),  axial  and  radial  velocities, 
respectively;  U in  is  the  inflow  velocity  in  the  propeller  fixed  frame. 


The  velocities  Vt  and  Va  are  computed  by  using  a  method  similar  to  that  of 
Greeley  and  Kerwin  [16],  where  a  tip  vortex  with  finite  core  is  utilized.  The  core 
radius  and  the  radial  contraction  angle  of  the  tip  vortex  at  the  blade  are  given 
from  experimental  information.  On  the  other  hand,  is  computed  by  using 
the  superposition  of  the  radial  velocities  induced  by  the  blade,  hub  and  wake 
singularities  resulting  from  applying  the  panel  method  on  the  conventional  grid. 


Figure  2-11:  The  geometry  of  grid  lines.  The  tip  vortex  detachment  point,  At^p,  is 
taken  downstream  of  the  actual  tip  due  to  the  contraction  of  the  wake. 


STEP  2 


@  Calculate  the  contraction  angle,  of  the  gridlines  along  the  blade  trailing  edge 
by  using  the  following  equation,  as  shown  in  Figure2-10. 


V 


M 

WAKE 


(2.3) 


with  and  Vr  computed  in  STEP  1. 


®  Find  the  location  of  the  tip  vortex  detachment  point  (also  called  the  “compu¬ 
tational”  tip),  Atip,  shown  in  Figure  2-11.  At  first,  the  tip  is  determined  by 
searching  among  the  streamlines  (corresponding  to  for  the  one  which 

starts  at  the  largest  blade  radial  location  and  does  not  intersect  the  propeller 
blade.  Due  to  the  contraction  of  the  wake,  the  location  of  the  computational 
tip  will  be  downstream  from  the  actual  tip.  The  details  of  the  location  of  the 
tip  will  be  discussed  in  the  following  chapters. 

®  Construct  the  flow  adapted  grid  on  the  blade.  Having  determined  the  location 
of  the  tip  vortex  detachment  point,  Aup,  the  grid  on  the  blade  is  determined, 
by  modifying  the  algorithm  on  which  the  blade  orthogonal  grid  was  based,  as 
follows: 

First,  the  arclengths  from  the  hub  to  the  tip  along  the  blade  leading  edge  and 
trailing  edge  are  divided  into  M  half-cosine  intervals,  with  M  being  the  number 
of  “spanwise”  panels  as  shown  in  Figure  2-11.  The  corresponding  arclengths 
are  given  as 
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Atip')  COS  /?n 


for  m  =  1,  2, 
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where  aie{m),  afe(m)  are  arc  length  from  the  leading  edge  of  the  hub  to  Aiim) 
and  from  the  trailing  edge  of  the  hub  to  Arim),  respectively,  with  defined 


Figure  2-12:  FLow  Adapted  Grid  on  N4990  propeller  blade;A’  =  40,  M  =  20. J  = 
1.270. 

In  the  case  of  the  blade  orthogonal  grid,  each  Ai{m)  is  connected  with  Arim) 
via  B-spline  curves  with  four  vertices  which  are  normal  to  the  blade  outline. 
In  the  present  case,  the  grid  lines  are  still  normal  to  the  blade  outline  at  the 
leading  edge,  but  now  they  form  an  angle  at  the  trailing  edge,  which  is  equal  to 
the  corresponding  wake  contraction  angle  9,  defined  earlier.  The  arcs  of  these 


grid  lines  are  then  divided  into  N/2  full-cosine  spaced  intervals,  with  N  being 
the  number  of  “chordwise”  panels. 


Figure  2-13:  Effect  of  wake  geometry  on  circulation  predicted  by  BEM  for  three 
hydrofoils  with  elliptic  chord  distribution  along  the  span;  [r/cj^ox  =  0.2,  a  =  5.73°, 
aspect  ratio  AR=3.  45°  backward  sweep  (top)  no  sweep  (middle)  45°  forward  sweep 
(bottom).  The  corresponding  flow  adapted  grids  are  also  shown.  Contraction  of  wake 
is  approximated  with  that  due  to  wing  thickness  effect. 


The  grid  on  both  sides  of  the  propeller  blade  is  then  determined  by  moving  the 
grid  of  the  blade  normal  to  the  three  dimensional  blade  camber  surface  by  an 
amount  equal  to  ±t/2. 
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•  Construct  the  new  trailing  wake  grid  based  on  the  contraction  angle,  6,  at  the 
tip.  The  ultimate  wake  geometry  is  kept  the  same  to  that  of  the  conventional 
grid. 
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Figure  2-14:  Circulation  distribution  on  a  circular  planform  hydrofoil;  [r/cj^ai  =  0.2, 
a  =  5.73°.  Predicted  by  applying  the  BEM  on  the  flow  adapted  grid;  before  and  after 
applying  the  IPK  condition. 


STEP  3 

•  Apply  the  panel  method  on  the  FLAG  which  was  generated  in  STEP  2.  De¬ 
termine  the  perturbation  potential  distribution  and  calculate  the  contraction 
angle  at  the  trailing  edge. 


Repeat  STEPS  2  and  3  until  the  contraction  angle  is  converged.  Usually  2  or  3 
iterations  are  enough.  Finally,  calculate  spanwise  circulation  distribution,  pressure 
coefficient  and  thrust/torque  on  the  blade.  This  grid  for  propeller  N4990  is  shown 
in  Figure  2-12.  This  iterative  procedure  for  determining  FLAG  has  been  applied  for 
several  wings  and  propellers  and  the  results  will  be  shown  in  the  next  section. 
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Figure  2-15:  The  effects  of  thickness  sources/or  sinks  on  the  flow  field  in  the  wakes 
of  elliptic  planform  hydrofoils.  The  effect  of  sinks  on  the  wake  is  stronger  than  that 
of  sources. 


2.4.2  The  FLAG  on  Wings 

The  BEM  is  applied  on  three  hydrofoils  with  the  following  characteristics: 

®  Elliptic  chord  distribution  in  the  spanwise  direction. 

®  Aspect  Ratio  AR=3. 
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•  Modified  NACA66  thickness  distribution  with  20%  maximum  thickness  to  chord 
ratio  (constant  in  the  spanwise  direction) 

•  Sweep  angle  varying  from  -45°  (forward  sweep)  to  0°  and  45°. 


s  /s 
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Figure  2-16:  Tip  vortex  trajectories  predicted  from  analysis  and  measured  in  experi¬ 
ment.  Elliptic  planform  hydrofoil;  a  =  15.5°,  [r /cj^ai  =  0.15,  aspect  ratio  AR  =  3. 

The  circulation  distributions  predicted  from  applying  the  BEM  on  the  conven¬ 
tional  and  the  flow  adapted  grid  are  shown  in  Figure  2-13.  Notice  that  as  the  sweep 
angle  decreases  the  difference  between  the  two  circulation  distributions  increases. 
This  is  due  to  the  drastic  contraction  of  the  wake  geometry,  also  shown  in  Figure 
2-13,  as  the  planform  is  swept  forward.  The  wake  contraction  is  increasing  as  the 
planform  is  swept  forward,  because  the  thickness  sinks  at  the  aft  part  of  the  plan- 
form,  shown  in  Figure  2-15,  have  a  stronger  effect  on  the  wake  streamlines  at  the 
tip  than  the  sources  at  the  forward  part.  In  addition,  the  circulation  distributions 
before  and  after  applying  the  IPK  condition  for  a  circular  wing  planform  hydrofoil 
are  shown  in  Figure  2-14. 


o - VLM  WITH  THICKNESS  CORRECTION 


Figure  2-17:  Circulation  distributions  predicted  from  BEM  (after  IPK  condition)  and 
Vortex-lattice  Method  (including  the  thickness  loading  coupling)  for  different  grid 
arrangements.  The  analytical  solution  of  zero  thickness  (Jordan,  73)  is  also  shown. 
Circular  planform  hydrofoil  ;  [r/c]mai  =  0.2,0;  =  5.73°.  Conventional  grid(top), 
BOG(middle),  FLAG(bottom). 
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Notice  that  the  two  curves  are  closer  to  each  other  than  they  were  in  the  case  of 
the  conventional  and  the  blade  orthogonal  grid  as  shown  in  Figures  2-4,  2-8. 


Figure  2-18:  Pressure  coefficients  predicted  from  BEM  applied  on  the  flow  adapted 
grid  (after  IPK  condition)  ;  propeller  N4119,  J  =  0.833,  where  r/R  is  defined  at  the 
trailing  edge. 

Finally,  the  elliptic  planform  hydrofoil  used  in  the  experiment  by  Arndt  et  al.  [2] 
is  analyzed  by  the  present  method.  In  particular,  the  trajectory  of  the  tip  vortex 
is  determined,  by  considering  only  the  effect  due  to  hydrofoil  thickness.  The  span- 
wise  coordinates  of  this  trajectory  are  then  superimposed  to  those  of  the  tip  vortex 
trajectory  predicted  by  Krasny  [37],  who  only  included  the  effects  of  the  wake  sheet 
roll-up.  The  results  are  shown,  together  with  the  experimental  results,  in  Figure  2- 


46 


16.  Notice  that,  under  the  examined  condition,  the  effects  of  the  wake  sheet  roll-up 
and  that  of  the  hydrofoil  thickness  on  the  predicted  tip  vortex  trajectory  are  equally 
important,  and  that  when  both  are  included,  the  trajectory  of  the  tip  vortex  matches 
the  observed  very  well. 


Figure  2-19:  Mean  velocity  vectors  along  the  trailing  edge  of  the  propeller  N4119  ; 
J  =  0.833.  Predicted  by  applying  the  BEM  on  the  flow  adapted  grid. 

However,  in  order  to  validate  the  prediction  methods  completely,  more  compar¬ 
isons  with  experiments  at  several  conditions  are  required.  In  addition,  further  com¬ 
putations  which  also  include  the  effects  of  viscosity  on  the  hydrofoil  loading,  thus  on 
the  strength  of  the  wake  vortex  sheet,  must  be  carried  out. 
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When  the  flow  adapted  grid  was  applied  on  a  circular  planform  hydrofoil,  the 
performance  of  the  boundary  element  method  was  found  to  improve  substantially.  In 
addition,  the  results  from  applying  a  vortex-lattice  method  on  the  same  grid  were 
found  to  be  in  very  good  agreement  to  those  from  the  boundary  element  method  as 
shown  in  Figure  2-17  [36]. 


Figure  2-20:  Circulation  distribution  on  the  propeller  N4119;  J  =  0.833.  Predicted 
by  applying  the  BEM  on  the  flow  adapted  grid;  before  and  after  applying  the  Iterative 
Pressure  Kutta  condition. 

Thus,  previous  differences  between  the  results  from  the  boundary  element  and  the 
vortex-lattice  method  have  been  reconciled  when  the  flow  adapted  grid  was  incorpo¬ 
rated.  In  this  comparison,  the  Vortex-Lattice  Method  (VLM)  is  modified  to  include 
FLAG  as  well  as  the  effect  of  coupling  between  thickness  and  loading  from  Kinnas 


[34]. 


2.4.3  The  FLAG  on  Propellers 

The  BEM  with  FLAG  is  applied  to  several  propeller  blades  and  the  results  are  com¬ 
pared  with  experimental  data. 


Figure  2-21:  Circulation  distributions  predicted  by  the  BEM  applied  on  the  conven¬ 
tional  and  the  flow  adapted  grid. 


®  PROPELLER  N4119 

N4119  is  a  typical  propeller  with  large  tip  chord  and  no  skew.  The  resulting 
pressure  distributions  are  shown  in  Figure  2-18.  In  this  figure  it  can  be  seen  that 
the  pressure  distribution  near  the  tip  region  obtained  by  using  the  flow  adapted 
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grid  is  less  singular  than  that  from  the  conventional  grid,  shown  in  Figure  2-5. 
The  resulting  total  mean  velocity  vectors  along  the  trailing  edge  are  shown  in 
Figure  2-19.  Comparing  Figure  2-19  to  2-6,  we  notice  that  the  singular  behavior 
of  the  velocity  vectors  at  the  tip  has  disappeared  and  that  the  resulting  flow  in 
the  wake  is  aligned  with  the  wake  gridlines.  Figure  2-20  shows  the  predicted 
circulation  distributions  before  and  after  the  iterative  pressure  Kutta  condition 
to  be  practically  identical.  In  the  case  of  the  conventional  grid,  the  difference 
between  the  circulation  distributions  (before  and  after  IPK),  shown  in  Figure 
2-3,  appears  to  be  slightly  larger.  In  addition,  the  number  of  iterations  for  the 
iterative  pressure  Kutta  condition  to  converge  (about  4-5)  is  much  smaller  than 
that  in  the  case  of  conventional  grid.  The  circulation  distributions  from  the 
conventional  and  the  flow  adapted  grids  are  shown  in  Figure  2-21.  They  appear 
to  be  very  different,  not  only  at  the  tip  but  over  a  large  part  of  the  propeller 
radius. 


Number  of  Blades  :  5 

Hub/Diameter  Ratio  :  0.3 

Section  Meanline  :  NACA  a=0.8 

Section  Thickness  Form  :  NACA66  (Modified) 


Table  2.1:  The  geometry  of  the  propeller  N4990. 
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Notice  that  the  values  of  the  circulation  distribution  from  FLAG  extend  up  to 
the  radius  of  the  computational  tip,  which  in  this  case  was  at  0.95i2.  Had  the 
circulation  distribution  from  FLAG  been  scaled  with  respect  to  its  correspond¬ 
ing  tip  radius,  it  would  result  into  a  circulation  distribution  which  would  be 
practically  identical  to  that  from  the  conventional  grid. 


Figure  2-22:  Panel  arrangement  for  propeller  N4990. 

®  PROPELLER  N4990 

N4990  is  a  highly  skewed  propeller.  The  geometry  of  propeller  N4990  is  given  in 
Table  2.1.  The  panel  arrangement  on  the  blades  and  hub  are  shown  in  Figure  2- 
22.  This  propeller  combines  a  wide  tip  geometry  and  a  high  skew  at  the  tip.  The 
pressure  distribution,  predicted  from  the  BEM  applied  on  the  conventional  and 
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the  flow  adapted  grids,  are  shown  in  Figures  2-23  and  2-24,  respectively.  Notice 
the  singular  behavior  of  the  pressure  at  the  tip  in  the  case  of  the  conventional 
grid.  It  should  be  noted  that  the  circulation  distributions  (not  shown)  from  the 
two  grids  are  almost  identical. 


Figure  2-23:  Pressure  coefficients  predicted  from  BEM  applied  on  the  conventional 
grid  (after  IPK  condition);  propeller  N4990,  J  —  1.270. 

This  is  attributed  to  the  negligible  change  in  the  radial  location  of  the  com¬ 
putational  tip,  due  to  the  combination  of  high  skew  (which  may  be  seen  as 
equivalent  to  the  “backward”  sweep  in  the  case  of  wings  as  shown  by  Kinnas 
et  al.  [36]  )  and  wide  chord  at  the  tip.  The  convergence  of  inviscid  values  of 
Kt,  Kq  and  r)  with  number  of  panels  is  shown  in  Table  2.2  for  the  design  J. 
It  should  be  noted  that  the  effect  of  hub  was  not  included  in  this  convergence 


test.  A  60  X  30  grid  seems  to  be  adequate  for  this  high  skew  geometry  case.  In 
order  to  compare  the  results  to  those  from  experiment  (supplied  by  Dr.  Jessup 
of  DTMB)  we  carried  out  a  calculation  with  a  60  x  30  grid  in  which  the  paneling 
on  the  hub  was  also  included. 


Figure  2-24:  Pressure  coefficients  predicted  from  BEM  applied  on  the  flow  adapted 
grid  (after  IPK  condition);  propeller  N4990,  J  =  1.270. 


The  viscous  eflPects  on  the  forces  were  estimated  via  a  constant  friction  coefficient 
Cp,  assumed  to  be  equal  to  0.005.  The  viscous  force  was  computed  by  the 
following  equation  ; 


Fv  —  -p  Cp  [  \Vtotal  \  V  tot  Aids 
2t  J  S 
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where  Vtotal  is  the  total  velocity  on  the  surface  of  the  blade.  Table  2.3  shows 
Kt,Kq  and  77  with  and  without  viscous  effects,  versus  the  measured  values. 
The  differences  between  the  experimental  and  the  computational  results  are 
within  5%. 


PROPELLER  N4990  (J=1.270) 

(  Cp  =  0.0&:  With  Leading  Edge  Suction  ) 


CONVENTIONAL 

Kt 

Kq  X  10 

V 

40c  X  20s 

0.2319 

0.5519" 

0.8492 

60c  X  30s 

0.2406 

0.5759 

0.8442 

80c  X  40s 

0.2384 

0.5718 

0.8430 

90c  X  40s 

0.2379 

0.5693 

0.8448 

Table  2.2:  Convergence  of  Kp  and  Kq  (C^  =  0). 


PROPELLER  N4990  (J=1.270) 

(  With  Hub) 


CONVENTIONAL 

Kt 

Kq  X  10 

n 

Cp  =  0.0000 

0.2392 

0.5753 

0.8403 

Cp  =  0.0050 

0.2292 

0.6732 

0.6851 

EXPERIMENT 

0.243 

0.691 

0.709 

Table  2.3:  Kp  and  Kq  with  and  without  viscous  effect. 


GRID 

J 

Kt 

ERROR 

X 

0 

ERROR 

CONVENTIONAL 

0.833 

0.149 

2% 

0.306 

9% 

1.100 

0.033 

3% 

0.114 

8% 

FLAG 

0.833 

0.143 

2% 

0.286 

2% 

1.100 

0.034 

0% 

0.113 

7% 

EXPERIMENT 

0.833 

0.146 

0.280 

1.100 

0.034 

0.106 

Table  2.4:  Force  coefficients  of  propeller  N4119  (  Cp  =  0.005  ). 
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These  force  parameters  were  compared  to  those  from  the  flow  adapted  grid  for 
different  J’s,  as  shown  in  Figure  2-25  and  Table  2.4.  The  coefficients  from  the 
flow  adapted  grid  were  much  closer  to  the  experimental  values.  Especially,  the 
value  of  Kq  was  much  closer  to  the  experimental  value  as  shown  in  Table  2.4. 


K. 


0.20 


0.15 


0.10 


0.05V 


0.00 


—EXPERIMENT 

-^ORIGINAL 

-^FLAG 

J  =  0.833 


J.  =  1.1 


4-.-.  ■  i  1. 


J 


0.6  0.8  1.0  1.2  1.4  s 


Figure  2-25:  ;  Force  coefficients  of  propeller  N4119.:  with  hub  effect,  Cp  =  0.005, 
J  =  0.833, 1.1. 

2.4.4  Numerical  Validation 

In  this  section  in  order  to  validate  the  results  of  panel  method  with  FLAG,  a  conver¬ 
gence  test  of  the  results  from  the  method  for  a  wide  tip  propeller  (N4119)  and  the 
consistency  test  of  the  method,  also  for  a  wide  tip  propeller  (N4118)  are  performed. 

®  Convergence  Test  for  Propeller  N4119. 

In  this  test,  the  convergence  of  the  circulation  distributions  is  examined,  as  the 
number  of  panels  in  both  chordwise  and  spanwise  directions  are  increased.  The 
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convergence  of  the  resulting  circulation  distribution  from  applying  the  BEM 
on  FLAG  with  varying  numbers  of  chordwise  and  spanwise  panels  is  shown  in 
Figure  2-26.  A  40  x  20  grid  appears  to  produce  convergent  results. 


Figure  2-26:  Convergence  test  for  N4119  propeller. 

•  Consistency  Test  for  Propeller  N4118 

A  test  of  the  consistency  of  lifting  surface  and  panel  method  can  be  made  by 
linearly  extrapolating  to  zero  thickness  the  results  of  a  sequence  of  panel  calcu¬ 
lations  for  different  thickness  ratios.  The  panel  calculations  are  made  for  a  set 
of  geometries  which  are  identical  except  for  a  constant  scale  factor  applied  to 
the  thickness.  The  circulation  distribution  for  zero  thickness  is  then  determined 
by  linearly  extrapolating  with  thickness  the  circulation  from  the  panel  method. 
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This  zero  thickness  “panel  method”  result  is  then  compared  to  the  result  from 
the  lifting  surface  analysis  method  by  Greeley  and  Kerwin  [16].  The  results 
for  the  N4118  propeller  are  shown  in  Figure  2-27.  The  circulation  distributions 
from  the  panel  method  (applied  on  FLAG)  for  different  thickness  scale  factors 
(100%  corresponds  to  the  original  thickness  distribution)  appear  to  extrapo¬ 
late  smoothly  (linearly)  to  the  circulation  distribution  from  the  lifting  surface 
method. 


Figure  2-27:  Consistency  test  for  N4118  propeller. 

Notice  that  as  the  thickness  decreases  the  location  of  the  computational  tip 
approaches  the  actual  tip.  The  extrapolated  panel  method  result  versus  the 
lifting  surface  result  is  shown  at  the  bottom  part  of  Figure  2-28.  The  consistency 
test  in  the  case  of  applying  the  panel  method  on  the  conventional  grid  has  been 
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performed  by  Hsin  [24]  and  the  results  are  shown  at  the  top  part  of  Figure  2-28. 
A  3%  difference  between  the  maximum  circulation  from  the  panel  method  and 
the  lifting  surface  method  was  reported  by  Hsin  [24] .  This  difference  appears 
to  become  smaller  in  the  case  of  the  flow  adapted  grid,  as  can  be  seen  in  Figure 
2-28. 

In  this  chapter,  the  flow  adapted  grid  was  incorporated  in  an  existing  BEM 
for  the  analysis  of  propeller  tip  flows.  The  proposed  grid  has  been  found  to 
predict  the  flow  and  pressures  at  the  tip  more  reliably  than  when  employing  the 
conventional  grid.  For  wide  tip  blades  with  zero  skew  (N4119)  the  predicted 
circulation  distributions  were  found  to  be  different  from  those  employing  the 
conventional  grid,  primarily  due  to  the  change  of  location  of  the  tip  detach¬ 
ment  point  in  the  case  of  the  flow  adapted  grid.  On  the  contrary,  for  wide  tip 
blades  with  high  skew  (N4990)  the  proposed  grid  did  not  affect  the  predicted 
circulation  distribution  even  though  in,  improved  the  predicted  pressures  at  the 
tip  substantially.  In  addition,  the  flow  adapted  grid  was  found  to  improve  the 
validity  of  the  consistency  test  between  the  boundary  element  and  vortex  lattice 
method  for  the  N4118  propeller. 


58 


Chapter  3 


The  Vortex  Sheet  Roll-Up 

In  the  previous  chapter,  the  trailing  wake  sheet  was  assumed  to  be  aligned  with  the 
velocities  at  the  trailing  edge  in  the  wake.  The  effect  of  wake  sheet  roll-up  was  ignored. 
In  reality,  the  trailing  wake  sheet  is  rolling  up  into  a  core  containing  the  vorticity  shed 
from  the  tip  region.  This  roll-up  effect  is  considered  to  be  of  the  primary  mechanism 
of  the  tip  vortex  structure.  In  addition,  when  the  propeller  is  heavily  loaded  or  the 
circumferential  velocity  due  to  the  rotation  of  the  propeller  is  large  compared  to  the 
forward  velocity  of  the  ship,  the  sheet  passes  very  close  to  the  following  blades  and 
the  wake  sheet  roll-up  results  in  a  change  in  loading,  velocity  field  and  pressures  on 
the  propeller  blades.  In  this  chapter,  the  vortex  sheet  roll-up  will  be  calculated  in 
two  and  three  dimensions  by  using  various  methods,  and  the  results  will  be  compared 
with  each  other.  In  Chapter  5,  the  wake  sheet  roll-up  will  be  included  in  the  flow 
adapted  grid. 

3.1  Two-Dimensional  Method 

3.1.1  Vortex  Blob  Method 

The  conjugate  complex  velocity  q{z)  induced  by  a  two-dimensional  vortex  sheet  of 
strength  7(s,  t)  =  dV/ds  situated  on  the  curve  C{s)  is  given  by  the  Rott-Birkhoff 
[67],  [6]  nonlinear  integro-differential  equation  ; 

1 
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(3.1) 


9(z) 


1  r  'y{s',t)ds' 
2z7r  Jc  z  —  z{s  it) 


+  Ue-  iVe 


where  Ug  and  Vg  are  the  components  of  the  incoming  flow  evaluated  at  z  as  shown  in 
Figure  3-1. 


(  U,V  ) 


Figure  3-1:  Velocity  and  vorticity  diagram. 

The  Cauchy  principal  value  is  assumed  for  the  integral  in  equation  (3.1)  in  order 
to  calculate  the  velocity  at  points  on  the  sheet.  Equation  (3.1)  has  been  generalized 
to  the  case  of  a  vortex  sheet  with  small  thickness  by  Moore  [57].  If  the  circulation  F 
is  chosen  as  the  Lagrangian  variable,  equation  (3.1)  may  be  written  from  Birkhoff  [6] 
as, 


dz 

dt 


(r.i) 


_L  /  dr' 

2i7r  Jc  z{Ti  t)  —  2:(r',  t) 


+  Ue-  iVe 


(3.2) 
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Equations  (3.1)  and  (3.2)  ensure  the  continuity  of  pressure  across  the  sheet  and  the 
conservation  of  circulation  around  a  segment  lying  between  any  two  points  moving 
with  the  sheet.  Equation  (3.2)  must  be  solved  numerically,  since  an  analytic  solution 
does  not  exist. 

In  the  vortex  blob  method,  the  vortex  sheet  is  replaced  by  a  finite  number  of 
points,  namely  “blobs” ,  which  are  assumed  to  be  rigid.  Then  the  motion  of  the  sheet 
is  approximated  by  calculating  the  trajectories  of  the  blob.  For  example,  when  the 
trailing  vortex  sheet  of  a  wing  is  represented  by  an  array  of  N  line  vortices,  equation 
(3.2)  reduces  to  an  initial  value  problem,  consisting  of  a  set  of  2N  first  order  ordinary 
differential  equations  whose  solution  requires  suitable  smoothing  technique.  Many 
numerical  attempts  have  been  made  with  different  smoothing  techniques,  such  as 
amalgamation  by  Ham  [17]  and  Fink  and  Soh  [13],  subvortex  model  by  Maskew  [46], 
rediscretization  by  Fink  and  Soh  [14]  and  Sarpkaya  [69],  and  cut-off  scheme  by  Krasny 
[37]  and  Kuwahara  [39]. 

In  this  section,  Krasny’s  method  [37]  will  be  described  briefly.  The  discretized 
form  of  the  equation  (3.2)  may  be  expressed  as. 


2^7r  ^  Zk-  Zj 


(3.3) 


In  desingularizing  the  vortex  sheet  equation  (3.3),  an  artificial  smoothing  param¬ 
eter,  S,  is  applied  at  the  right-hand  side  of  equation  (3.3). 


Uk-lVh  =  7^2^  - 


Zk-Zjf 


2myzk-Zj  \zk-Zj\ 


+  (C/e)fc  -  i{Ve)k 


(3.4) 


The  equation  (3.4)  is  no  longer  singular  at  Zk  =  Zj  for  5  >  0.  As  ^  goes  to  zero, 
equation  (3.4)  converges  to  the  original  equation  (3.3).  Krasny  [37]  demonstrated 
that  the  numerical  solutions  of  the  vortex  blob  equations  (3.4)  for  different  values  of 
d  converge  to  a  limit  curve  as  5  0  and  this  curve  can  be  interpreted  as  a  weak 


solution  of  the  original  vortex  sheet  evolution  equation  (3.1). 


Figure  3-2:  Roll-up  behind  a  lifting  line  with  elliptic  loading  predicted  by 
Krasny(1987)  ;  (5  =  0.05,  N  =  200,  At  =  0.01. 

Two  examples  are  shown  in  Figures  3-2  and  3-3.  The  cases  studied  are  the  vortex 
sheet  roll-up  for  a  lifting  line  with  elliptic  loading  and  an  embedded  vortex  sheet  in 
two-dimensional  uniform  flow.  The  methods  predict  the  vortex  sheet  roll-up  very 
smoothly  even  in  the  late  stage  because  5  provides  the  damping  of  short  waves.  The 
method  also  shows  good  convergence  with  number  of  panels.  However,  the  use  of 
vortex  blobs  has  the  following  problems  : 

®  Mathematically,  the  linear  sum  of  vortices  cannot  constitute  an  exact  solution 
of  the  nonlinear  equations  of  motion. 
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•  Near  the  core,  each  blob  is  not  allowed  to  distort  because  the  blob  is  assumed 
to  have  an  invariable  core  shape  and  size. 

•  Finally,  it  is  time  consuming  and  difficult  to  extend  to  three  dimensions. 


Figure  3-3:  A  vortex  sheet  embedded  in  two-dimensional  uniform  flow;  5  =  0.5,  N  = 
400,  At  =  0.1.  Predicted  by  Krasny  [37]. 


3.1.2  Discrete  Vortex  Method 

The  velocity  induced  by  the  vorticity  concentrated  in  a  bounded  region  is  given  by 
Bachelor  [5] 
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Figure  3-4:  Cross  section  of  the  wake  sheet  behind  a  lifting  line  of  2  ft  span  with 
elliptic  loading.  The  circulation  at  midspan  is  Sft^/sec.  99  trailers  are  taken  at  even 
spacing. 

If  the  vorticity  is  concentrated  in  a  single  curved  filament  of  circulation  F,  equation 
(3.5)  reduces  to  the  Biot-Savart  law 
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(3.6) 


r  r  {r{s)-r'{s))  X  dr'js) 

^  47r/c  \r{s)  -  r' {s')f 

where  r(s)  describes  a  position  vector  of  the  filament  centerline  in  terms  of  the  ar- 
clength  s. 

Equation  (3.6)  yields  an  infinite  self-induced  velocity  if  the  filament  is  curved  and 
zero  induced  velocity  if  it  is  straight  as  explained  by  Bachelor  [5].  Several  methods 
have  been  introduced  to  handle  this,  such  as  the  Local  Induction  Approximation  by 
Hama  [18].  However,  the  only  way  to  get  the  desired  shape  of  the  roll-up,  is  to  ignore 
this  term. 

In  this  section,  the  vortex  sheet  is  represented  by  a  series  of  discrete  line  vortices. 
The  calculations  are  done  midway  between  vortices  to  avoid  an  infinite  self-induced 
velocity.  For  each  location,  the  velocity  due  to  each  element  is  calculated  and  summed 
to  obtain  the  total  induced  velocity  V  at  each  point  on  the  sheet.  Starting  at  the 
trailing  edge,  a  new  location  of  the  sheet  is  formed  by 

®(n,p)  ~  ®(n,p— 1)  T  ^(n,p— 1)  ‘  At. 

where  n  describes  the  number  of  iteration  and  p  is  the  index  for  station  along  the 
streamwise  direction.  The  result  for  the  lifting  line  with  elliptic  loading  is  shown  in 
Figure  3-4.  This  method  seems  to  predict  the  shape  of  the  roll-up  smoothly.  The 
method  can  very  easily  be  extended  to  three  dimensions.  However,  the  choice  of  vortex 
spacing  has  an  effect  on  the  ultimate  shape  of  the  roll-up  and  also,  there  are  difficulties 
with  this  method  even  for  two-dimensional  flow  simulation.  First,  the  vortex  filaments 
are  singularities  and  hence  create  large  velocities  in  their  neighborhood.  This  causes 
instabilities  and  physically  impossible  sheet  crossings  along  and  near  the  edges  of  the 
sheet.  The  second  difficulty  is  in  the  CPU  per  time  step.  The  number  of  operations 
required  for  the  velocity  calculation  is  proportional  to  N'^,  where  N  is  the  number 
of  vortices.  Thus,  the  CPU  time  increases  significantly  as  the  number  of  vortices 
increases. 
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3.1.3  Panel  Method 


There  have  been  a  large  number  of  attempts  to  model  the  vortex  sheet  motion  by 
replacing  the  continuous  vortex  sheet  with  a  finite  number  of  discrete  vortices.  How¬ 
ever,  it  appears  that  despite  these  various  attempts,  the  discrete  vortex  representation 
results  inevitably  in  numerical  instabilities,  i.e.  chaotic  motion,  when  the  number  of 
vortices  is  increased.  It  must  also  be  noticed  that  the  question  of  the  existence  and 
uniqueness  of  the  solutions  for  equations  describing  vortex  sheet  motion  has  not  yet 
been  fully  resolved  as  mentioned  by  Moore  [58].  So,  it  can  be  said  that  the  discrete 
vortex  method  is  not  adequate  to  compute  vortex  sheet  roll-up  reliably.  Recently, 
panel  methods  have  been  applied  in  order  to  get  more  reliable  and  accurate  rep¬ 
resentations  of  the  vortex  sheet  roll-up  in  two  dimensions.  The  method  has  been 
extensively  described  by  Faltinsen  and  Petterson  [12]  and  Hoeijmakers  [22].  The  ear¬ 
lier  attempt  by  Mokry  and  Rainbird  [55]  using  a  low-order  panel  method  seemed 
promising,  but  numerical  instabilities  still  appeared.  It  has  been  suggested  that  the 
use  of  a  high-order  panel  method  is  more  accurate  than  the  low-order  panel  method 
or  the  discrete  vortex  method  to  compute  the  velocity  field,  since  it  precludes  the 
appearance  of  instabilities  in  the  vortex  sheet  due  to  the  spurious  numerical  effects 
introduced  by  a  too  crude  representation. 

In  this  section,  a  high  order-panel  method  is  used  to  compute  the  velocity  field 
induced  by  the  multiple  segmented  vortex  sheet  in  two  dimensions.  The  vortex  sheet 
is  modeled  by  piecewise  continuous  dipole  distributions.  The  velocity,  q,  at  a  point 
{x,y)  induced  by  a  dipole  distribution  r(t)  on  a  sheet  C{t)  is  written  as 

where  t  is  a  parametric  variable  along  the  sheet  C  {t)  and  dl  is  a  differential  element 
of  length  along  C{t). 

The  curve  C[t)  is  divided  into  number  of  panels  as  shown  in  Figure  3-5.  Then, 
equation  (3.7)  becomes 
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(3.8) 


Qi  =  ^(f){xi,yi)  = 


—  f  dT{tj)  ^ 

27r  j  Jc(tj)  dt 


dlj 


where  subscript  i  is  an  index  for  a  field  point  and  j  denotes  the  functions  on  the  j-th. 
segment.  On  each  panel,  the  function  r(t)  is  approximated  by  piecewise  quadratic 
representations.  Let  r(t)  =  cio  +  o-it  + 


I 


Figure  3-5:  Vortex  sheet  arrangement. 


Then,  equation  (3.8)  can  be  expressed  as  : 


q,  =  V<^(a;i,y,)  =  -^^^^^^(ai-H2a2t,)v|^^logr,jj  dl^  (3.9) 

The  integrals  in  equation  (3.9)  can  be  expressed  in  closed  form.  For  a  field  point 
in  the  far  field,  the  computation  time  is  saved  by  using  the  multipole  expansion 
rather  than  the  closed  form.  To  avoid  spurious  effects  from  the  panel  edges  where  the 
geometry  and  the  dipole  distribution  may  be  discontinuous,  the  panel  midpoints  are 


chosen  as  the  control  points. 

The  detailed  computational  scheme  is  as  follows. 

1.  For  a  given  position  and  strength  of  the  sheet  at  time  r,  the  induced  velocity 
components  {Uj,  Vj)  are  computed  at  the  panel  midpoints  {tj). 


VORTEX  BLOB 
PANEL  METHOD 


Figure  3-6:  The  shape  of  the  roll-up  behind  a  lifting  line  with  elliptic  loading. 

2.  The  panel  midpoints  are  advanced  in  time  by  a  simple  Euler  scheme,  i.e.  by 
an  amount  of  T^Ar),  where  Ar  is  chosen  such  that  each  midpoint  is 

not  moved  more  than  Atj.  The  new  position  of  the  end  points  ty  and  are 
found  by  a  quadratic  extrapolation  from  the  new  values  of  the  three  nearest 
midpoints. 
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3.  The  new  values  of  x{ij),  y{ij)  and  r{ij)  are  then  the  input  data  for  a  cubic 
spline  procedure. 


Figure  3-7:  Convergence  test  of  the  results  from  the  panel  method.  The  vortex  sheet 
at  t  =  1.0  behind  the  lifting  line  with  elliptic  loading. 

4.  With  the  cubic  spline  representation  for  x{t),  y{t)  and  T{t),  a  rediscretization 
scheme,  called  “adapted  curvature- dependent”  [22],  is  applied,  which  computes 
the  values  of  the  panel  edges  (tj).  Using  the  rediscretization  method,  the  panel 
edge  points  are  computed  such  that  tj+i  =  tj  +  Amax  unless  Atj  spans  an  arc 
of  more  than  9max  degrees  of  a  circle  with  a  radius  equal  to  the  average  radius 
of  curvature  within  [tj,tj+i].  In  the  latter  case,  tj  =  9max/kn,  where  is  the 
average  radius  of  curvature  in  the  interval.  The  value  of  A^ax  and  9max  are  user 
specified.  For  example  in  this  chapter,  Amax  and  9max  are  set  equal  to  the  ratio 
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of  total  arclength  of  the  sheet  to  the  number  of  panels  and  45°,  respectively. 

5.  Finally,  the  new  values  of  y{tj)  and  r(tj)  are  computed  via  the  cubic 

spline  approximation. 


Figure  3-8;  The  wake  profile  for  a  heaving  hydrofoil. 

This  procedure  is  repeated  until  the  geometry  of  the  vortex  sheet  is  converged. 
As  a  test  of  this  method,  the  same  problem  as  that  shown  in  Figure  3-2  is  solved  and 
the  result  is  compared  to  that  from  applying  the  vortex  blob  method,  described  in 
Section  3.1.1.  The  resulting  geometries  of  the  vortex  sheet  roll-up  predicted  by  the 
two  methods  are  shown  in  Figure  3-6.  The  comparison  shows  that  the  results  from 
the  two  methods  agree  well  to  each  other  except  inside  the  core.  This  discrepancy 
in  the  core  will  be  discussed  in  section  3.2.3.  A  convergence  test  of  the  results  from 
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this  method  is  shown  in  Figure  3-7.  Figure  3-7  shows  a  clear  convergence  trend  with 
increasing  number  of  panels.  This  method  is  also  applied  to  the  two-dimensional 
unsteady  hydrofoil  problem  for  the  purpose  of  extension  to  the  three-dimensional 
unsteady  propeller  problem. 


Figure  3-9:  Comparison  of  analytic  and  numerical  solutions  for  the  heaving  hydrofoil. 

Figure  3-8  shows  the  wake  shape  behind  a  heaving  hydrofoil.  In  Figure  3-9,  it 
is  shown  that  the  lifting  coefficient  from  the  present  panel  method  agrees  very  well 
with  that  from  the  analytic  method,  given  by  the  well  known  “Theodorsen  function” 
(Newman  [62]).  For  a  two-dimensional  calculation  of  the  vortex  sheet  roll-up,  it  is 
shown  that  the  high-order  panel  method  predicts  the  geometry  of  the  vortex  sheet 
roll-up  reliably  and  smoothly.  In  addition,  this  method  is  much  faster  than  the  vortex 
blob  method  and  is  very  easy  to  extend  to  three  dimensions. 
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3.2  Three-Dimensional  Methods 


A  number  of  three-dimensional  flow  models  have  been  implemented  for  the  prediction 
of  wake  sheet  roll-up.  These  consist  of  the  discrete  vortex/strip  theory  methods, 
vortex  lattice  methods  and  panel  methods.  Existing  models  will  be  summarized  first, 
and  then  a  new  high-order  panel  method  will  be  described. 


X 


3.2.1  Discrete  Vortex/Strip  Theory  Method 

The  discrete  vortex  method  can  be  easily  extended  to  three  dimensions  in  a  strip- 
wise  sense.  This  method  can  predict  the  geometry  of  the  vortex  sheet.  However, 
the  method  does  not  include  three-dimensional  effects  such  as  the  concentration  of 
vorticity  in  the  core  region  because  the  streamwise  velocity  is  assumed  to  be  constant. 
This  three-dimensional  effect  was  first  modeled  by  Cummings  [8].  In  this  section,  his 
modeling  and  numerical  implementation  will  be  discussed.  The  classical  discrete  vor¬ 
tex  method  as  mentioned  earlier  suffers  from  lack  of  convergence.  The  ultimate  shape 
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computed  in  the  tip  region  depends  on  the  choice  of  the  vortex  spacing.  To  avoid 
this  problem,  a  model  of  the  tip  vortex  core  is  introduced.  The  tip  vortex  core  model 
assumes  that  the  tip  vorticity  is  spread  out  over  a  finite  thickness  (2a),  as  shown  in 
Figure  3-10,  where  a  is  a  boundary  layer  thickness  at  the  pressure  side  of  the  tip. 
Within  this  model,  even  if  the  spacing  becomes  fine,  the  resulting  shape  of  the  vortex 
sheet  will  not  be  singular  at  the  tip.  As  for  the  modeling  of  the  three-dimensional  ef¬ 
fect,  it  will  be  impossible  for  any  concentration  of  vorticity  (in  addition  to  that  due  to 
the  two  dimensional  roll-up)  to  occur  in  a  model  which  assumes  that  the  streamwise 
velocity  is  constant.  A  concentration  of  vorticity  implies  a  concentration  of  fluid  since 
the  vortices  move  with  the  fluid.  The  only  model  which  can  produce  a  concentration 
of  fluid  is  a  sink  in  the  core,  which  increases  the  axial  velocity  and  decreases  the  core 
pressure.  Far  downstream,  a  concentration  of  vorticity  can  occur  in  two  ways.  The 
roll-up  of  the  vortex  sheet  increases  the  circulation  surrounding  the  core  and  any  flow 
into  the  core  brings  in  circulation  from  the  sheet.  The  first  effect,  due  to  the  vortex 
sheet  roll-up,  is  small  in  any  case,  because  the  sheet  vorticity  is  small  compared  to 
that  in  the  core.  The  second  effect  is  a  major  source  of  increasing  core  circulation  due 
to  absorbing  of  the  trailing  vortex  sheet  into  the  core.  This  effect  can  be  modeled  as 
follows.  A  length  Ax  of  the  core  is  considered  with  a  pressure  force  F  acting  on  it  in 
the  axial  direction  as  shown  in  Figure  3-11.  The  radial  velocity  V  is  assumed  at  first 
to  be  zero.  Consider  the  continuity  and  momentum  equations  in  the  ^-direction  over 
a  control  volume  of  radius  a  and  length  Ax.  As  explained  in  Appendix  B,  a  can  be 
assumed  to  be  constant  and  the  flow  inside  the  core(as  defined  in  Appendix  B)  to  be 
inviscid. 

From  the  continuity  equation, 

2V 

U2  =  Ui  - -  Ax  (3.10) 

a 

where  Uq  is  an  inflow  velocity  outside  of  the  core;  ui  and  Fi  are  the  induced  velocity 
and  the  acting  force  on  the  upstream  side  of  the  core,  respectively;  U2  and  F2  are  the 
same  variables  on  the  downstream  side  of  the  core. 


From  the  momentum  equation  in  the  x-direction, 

pTxc?u\  +  Fi  —  pTxa^u\  +  27i:'apVUoAx  +  F2 

F  =  Fi  -  F2  =  pwa  (auj  -  aul  -  2AxVUo)  (3.11) 

Combining  equation  (3.10)  and  (3.11), 


a 

L  x2  4T 

1/  = 

4Ax 

Uo  —  2ui  -f  W(17o  '2ui)  -1 

V  pTT^a^ 

V 


Figure  3-11:  Control  volume  of  the  vortex  core. 

The  resulting  radial  velocity  V  means  that  a  sink  of  strength  27ray'  per  unit  length 
is  pulling  fluid  and  sheet  vorticity  with  it  toward  the  core.  In  this  calculation,  the 
stretching  effect  of  the  vortex  sheet  is  neglected.  This  model  predicts  a  decreased 
pressure  on  the  downstream  end  and  a  resultant  increase  of  axial  velocity  and  sink 
strength.  The  calculation  is  repeated  until  the  radial  velocity  converges  within  a 
desired  tolerance.  The  new  vortex  positions  may  now  be  calculated  and  the  process 
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moves  to  the  next  downstream  position.  Figure  3-12  shows  the  geometry  of  the  vortex 
sheet  roll-up  behind  an  elliptically  loaded  lifting  line  by  applying  the  core  and  sink 
model.  In  this  figure,  the  assumed  core  is  also  shown.  In  Figures  3-4  and  3-12,  it  is 
shown  that  the  concentration  of  the  vorticity  does  not  give  an  effect  on  the  geometry 
of  the  vortex  sheet.  Only  the  pressure  and  the  trajectory  of  the  tip  vortex  are  affected. 
The  calculation  of  pressure  in  the  core  is  explained  in  Appendix  B.  The  use  of  the 
strip  theory  model  in  the  discrete  vortex  method  provides  savings  on  the  computation 
time. 


Y 


Figure  3-12:  Cross  section  behind  an  elliptically  loaded  lifting  line  with  a  core  and 
sink  model.  Predicted  by  using  the  model  of  Cummings[8]. 


The  described  core  model  predicts  both  the  minimum  pressure  in  the  tip  vortex 
core  and  the  motion  of  the  sheet  with  reasonable  accuracy.  On  the  other  hand,  the 
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velocity  and  pressure  distributions  on  the  body  can  not  be  accurately  predicted  with 
this  model.  In  addition,  this  method  is  not  fully  three-dimensional. 

3.2.2  Vortex  Lattice  Method 

This  method  approximates  the  bound  vortex  sheet  by  a  bound  vortex  lattice  and  the 
free  vortex  sheet  by  a  set  of  segmented  free  vortex  lines.  Each  element  is  associated 
with  a  horseshoe  vortex.  Discrete  line  vortices  trail  to  infinite  downstream  along 
the  trailing  and  side  edges.  Their  alignment  with  the  local  flow  direction  provides 
the  necessary  boundary  conditions  on  the  vortex  sheet.  The  strengths  of  vortices 
are  determined  through  the  use  of  an  iterative  technique  by  imposing  the  kinematic 
condition  at  the  midpoint  of  the  3/4  chord  line  of  the  elements.  This  method  has 
been  widely  used  and  considerably  improved  by  Kandil  [30],  Rom  [65],  Almosnino  [1], 
and  Keenan  [31].  The  results  have  shown  that  the  model  is  rather  crude  and  exhibits 
undesirable  singular  behavior.  Even  though  overall  forces  are  predicted  well,  the  local 
velocity  or  pressure  distribution  on  the  body  and  the  vortex  sheet  geometry  are  not 
sufficiently  accurate.  In  addition,  increasing  the  number  of  vortices  makes  the  matter 
worse  as  shown  by  Rusak  et.  al  [68]. 

3.2.3  Panel  Method 

In  the  panel  method,  the  vortex  sheet  is  modeled  by  a  piecewise  continuous  dipole 
distribution,  thus  reducing  the  singular  behavior,  which  is  present  in  the  case  of  line 
vortices.  Many  attempts  have  been  made  with  this  method,  such  as  Suciu  and  Morino 
[72],  Johnson  et  al.  [29],  and  Hoeijmakers  [21].  In  this  section,  a  new  high-order  panel 
method  will  be  introduced. 

Assuming  that  the  circulation  is  given  as  r(t)  along  the  arclength  (t)  ,  then  the 
dipole  strength  on  the  vortex  sheet  can  be  found  from, 


7(t)  = 


+ 


—  u 
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Therefore 


^  ^  _  dA(j){t) 

dt  dt  dt 

r(t)  = 

where  f  is  a  parametric  variable  along  the  lifting  line,  and  u'^  and  u~  are  velocities 
on  the  upper  and  lower  surface,  respectively  as  shown  in  Figure  3-13. 


Figure  3-13:  Velocity  diagram  on  a  lifting  surface. 

Equation  (3.7)  can  be  generalized  to  the  three-dimensional  problem  by  replacing 
the  two-dimensional  source  expression  logr  and  2ii  by  the  three-dimensional  source 
expression  1/r  and  — 47r. 

,  =  V0(...)  =  i/^rwv(|i)  ds  (3.12) 

where  S  is  the  surface  of  the  vortex  sheet. 
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The  vortex  sheet  is  divided  into  a  number  of  panels  as  shown  in  Figure  3-14.  The 
control  point  is  set  at  the  centroid  of  the  panel.  Then,  equation  (3.12)  becomes 


dSj 


(3.13) 


INITIAL  GEOMETRY 


AFTER  2-ND  ITERATION 


AFTER  3-RD  ITERATION 


Figure  3-14:  The  geometry  of  the  trailing  vortex  sheet  for  a  lifting  line  with  elliptic 
loading;  Uq  =  1,  At  =  0.15  and  Fq  =  1.0/t^/sec.  The  low-order  panel  method  is 
used  with  16  and  31  spanwise  and  streamwise  number  of  panels,  respectively.  Notice 
divergence  of  roll-up  shape  at  the  3-rd  iteration. 
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If  a  planar  quadrilateral  panel  and  a  constant  strength  dipole  on  the  panel  are 
assumed,  equation  (3.13)  can  be  rewritten  as 


(3.14) 


The  integrals  in  equation  (3.14)  can  be  computed  analytically  by  using  the  ex¬ 
pressions  developed  by  Newman  [63]. 


Figure  3-15:  The  geometry  of  hyperboloidal  panel. 


The  velocities  at  the  control  points  on  the  vortex  sheet  are  computed  first  ,  and 
the  velocities  at  the  edges  of  the  panels  are  extrapolated  from  those  at  the  control 
points.  Then  the  panel  moves  with  velocity  {u  +  U,v  +  V^w  +  W),  where  {U,  V,  W) 
are  the  components  of  the  incoming  velocity. 


Figure  3-16;  Bi-quadratic  dipole  distribution  on  a  panel  with  nine  nodes. 

This  procedure  is  repeated  until  the  geometry  of  the  vortex  sheet  converges.  The 
result  for  a  lifting  line  with  elliptic  loading,  is  shown  in  Figure  3-14.  The  low-order 
panel  method  has  failed  after  three  iterations.  It  seems  that  the  planar  quadrilateral 
panel  can  not  model  the  highly  rolled-up  region  accurately,  due  to  the  discontinuity 
(gap)  in  geometry  and  in  strength  of  adjacent  panels.  As  a  more  accurate  and  reliable 
tool,  a  hyperboloidal  panel  and  a  high  order  strength  of  the  dipole  are  used. 

A  hyperboloidal  panel  is  defined  in  Figure  3-15.  Surface  E  is  a  hyperboloidal  sur¬ 
face  determined  by  four  vertices  ajij,  Xi^2,  ^2,1,  ^2,2,  and  the  surface  5  is  a  projected 
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nates 

{^,V)  =  ' 


The  vertices  thus  can  be  expressed  as  follows  ; 


^1,1 

1-1-1  1 

ao 

^1,2 

11-1-1 

ai 

^2,1 

1-11-1 

_  ^2.2  _ 

1111 

_  ^3  _ 

Figure  3-18:  Second  or  higher  iteration  of  the  scheme. 

By  inverting  this  matrix,  oq,  Ui,  02  and  03  can  be  obtained  from  the  vertices,  (Hsin 
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[24]). 


1-1-1  1 
1  1-1-1 
1-1  1-1 
1111 

A  bi-quadratic  dipole  strength  is  assumed  on  the  panel.  The  strength  of  the  dipole 
in  the  local  coordinate  system  can  then  be  expressed  as 

v)  —  bo  +  bi^  +  6277  -t-  bz^Tj  -f  641^^  -f  6577^  -l-  bo^tf  +  bj^^T]  -i-  (3.16) 

The  coefficients  bi  can  be  determined  from  the  dipole  strengths  at  the  nine  node 
points  on  the  panel,  as  shown  in  Figure  3-16.  The  9x9  matrix  equation,  which  must 
be  solved,  is 

1-11-1  1-11-11 
1  -1  1  0  0  0  0  0  0 

1-11  1-1  11-11 
1  0  0-1  0-11  0  0 

1  00  0  0  -1  0  00 

1  0  0  1  0-11  0  0 

1  11-1-1-11  11 
1  1  1  0  0  0  0  0  0 

1  11  1  1  11  11 

Now,  from  equations  (3.15)  and  (3.16),  equation  (3.13)  can  be  rewritten  as, 

V(f)i  =  —  XI  /  (^0  ^2^7  +  bz^rj  +  64^^  +  b^rf  +  be^rf  +  by^'^r] 
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The  calculation  of  the  integrals  in  equation  (3.17)  over  the  hyperboloidal  panel 
(Sj),  is  described  in  Appendix  C.  The  following  iteration  scheme  is  introduced  in  order 
to  satisfy  the  force- free  condition  on  each  panel  in  the  vortex  sheet.  To  accelerate 
the  iteration  process,  an  idea  similar  to  that  in  over-relaxation  schemes  is  adopted  in 
the  first  iteration  (Nagati  et  al.  [61]).  Namely,  in  the  first  iteration  the  calculation  of 
the  wake  surface  is  proceeding  from  the  first  station  to  the  last  section  downstream 
of  the  trailing  edge  as  follows  (shown  also  schematically  in  Figure  3-17). 


Figure  3-19:  The  geometry  of  vortex  sheet  behind  a  lifting  line  with  elliptic  loading. 
Predicted  by  the  present  method. 

0  The  total  velocities  at  the  control  points  of  the  current  row  of  the  panels  are 
computed  from  equation  (3.17).  The  velocity  at  the  edge  of  the  vortex  sheet  is 
computed  from  extrapolation. 
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•  Use  these  velocities  to  find  new  location  of  panels  on  the  current  row. 

•  The  panels  downstream  of  the  current  row  are  moved  by  the  same  amount  as 
those  at  the  current  row. 

•  Repeat  the  computation  until  the  last  row  of  panels. 

At  the  second  or  higher  iteration,  the  procedure  is  altered  due  to  the  fact  that 
the  geometry  of  the  upstream  vortex  sheet  will  be  aflFected  by  the  relocation  of  the 
vortex  sheet  downstream,  which  is  ignored  in  the  first  iteration.  The  following  steps 
are  used,  as  also  shown  in  Figure  3-18, 


Y  AT  Z=  1.0  X  SPAN 
At  =  n  m 


Figure  3-20:  Shape  of  the  vortex  sheet  with  number  of  iterations. 


•  Calculate  the  total  velocities  at  all  the  control  points  in  the  vortex  sheet. 


®  Move  all  the  panels  simultaneously  and  find  the  new  location  of  the  vortex 
sheet. 

®  Repeat  the  previous  steps  until  the  vortex  sheet  geometry  is  converged. 

To  maintain  the  accuracy  of  the  spatial  discretization  for  a  stretching  vortex  sheet, 
a  rediscretization  scheme  is  applied  as  the  roll-up  region  increases.  The  discretization 
scheme  has  to  be  refined  in  regions  where  large  variations  in  geometry  or  dipole 
distribution  appear  during  the  computation. 
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Figure  3-21:  Convergence  test  of  a  lifting  line  with  elliptic  loading  with  respect  to 
the  number  of  spanwise  panels. 


To  accommodate  these  requirements,  the  adapted  curvature-dependent  panel  scheme 
is  used  at  each  iteration,  similar  to  that  used  in  the  two-dimensional  panel  method  as 
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explained  in  section  3.1.3.  The  basic  panel  size  I^rnax  is  chosen  such  that  the  desired 
degree  of  accuracy  in  regions  where  no  large  variations  occur  is  ensured.  The  second 
parameter  dmax  then  ensures  that  at  highly  curved  parts  of  the  sheet  the  panel  size 
is  reduced  so  that  in  each  region  the  accuracy  is  maintained.  Usually  6max  is  chosen 
such  that  a  circle  having  a  radius  equal  to  the  radius  of  curvature  is  represented  by 
12-18  panels.  This  numerical  scheme  is  applied  to  a  lifting  line  with  elliptic  loading. 
The  resulting  geometry  of  the  vortex  sheet  roll-up  is  shown  in  Figure  3-19.  Note  that 
this  method  predicts  the  vortex  sheet  roll-up  very  smoothly  with  more  turns  (i.e.  at 
larger  distance  from  the  trailing  edge)  than  that  in  the  discrete  vortex  method. 


Y  AT  Z=  1.0  X  SPAN 
At  =n  ^0 


Figure  3-22:  Convergence  test  of  a  lifting  line  with  elliptic  loading  with  respect  to 
the  number  of  streamwise  panels. 

The  difference  among  the  iterations  is  shown  in  Figure  3-20.  This  figure  shows  that 


two  iterations  is  enough  usually  to  get  a  converged  result.  Finally,  more  convergence 
tests  are  shown  in  Figures  3-21  and  Figure  3-22.  Note  that  as  the  number  of  panels 
along  the  spanwise  and  streamwise  direction  increases,  the  geometries  of  the  vortex 
sheet  appear  to  converge  to  a  limit.  It  can  be  concluded  that  the  present  high-order 
panel  method  predicts  a  smooth  and  reliable  (convergent)  geometry  of  the  vortex 
sheet  roll-up  in  a  computationally  efficient  manner. 


Chapter  4 


The  FLAG  with  the  Wake  Sheet 
Roll-Up 


In  Chapter  2,  the  flow  adapted  grid  was  incorporated  in  the  panel  method  for  lifting 
surface  flows.  It  was  shown  that  using  the  flow  adapted  grid  in  the  panel  method, 
the  numerical  results  for  a  highly  skewed  propeller  and  a  propeller  with  a  large  tip 
chord,  especially  in  the  tip  region,  were  improved  substantially.  However,  the  flow 
adapted  grid  does  not  include  the  effect  of  a  wake  sheet  roll-up,  which  changes  the 
trajectory  of  the  tip  vortex  as  well  as  the  pitch  of  the  wake  sheet.  In  this  chapter,  the 
flow  adapted  grid  with  the  effect  of  the  wake  sheet  roll-up  will  be  incorporated  in  the 
panel  method.  In  calculating  the  wake  sheet  roll-up,  the  higher  order  panel  described 
in  the  previous  chapter  will  be  used.  The  flow  adapted  grid  will  be  constructed  by 
using  the  geometry  of  the  wake  sheet  roll-up,  in  an  iterative  sense. 

The  essential  elements  of  the  present  flow  model,  as  mentioned  in  Appendix  B,  are 
:  the  lifting  body,  the  trailing  wake  sheet,  the  sheet  emerging  from  the  tip  (sheath) 
and  the  roll-up  core  fed  by  the  tip  vortex  sheet.  On  each  of  these  elements,  the 
following  boundary  conditions  must  be  imposed: 

•  The  body  surface  is  impermeable. 

•  The  wake  sheet  cannot  support  a  pressure  difference  and  is  impermeable  as  well. 

•  The  Kutta  condition  is  imposed  along  the  trailing  edge  of  the  blade. 


90 


The  panel  method  will  be  modified  to  satisfy  the  boundary  conditions  on  the  body 
as  well  as  in  the  wake  sheet.  The  enhanced  panel  method  is  expected  not  only  to 
improve  the  accuracy  of  the  predicted  pressure  distributions  at  the  blade  tip,  but  also 
to  provide  the  foundation  for  predicting  the  tip  vortex  evolution. 

4.1  Mathematical  Formulation 


Figure  4-1:  A  control  volume  for  a  multidomain  problem. 

The  flow  is  assumed  to  be  inviscid,  incompressible  and  irrotational  everywhere,  except 
in  the  thin  wake  sheet.  The  flow  velocity  can  be  written  as  a  superposition  of  the 
incoming  flow  Uoo,  the  velocities  induced  by  the  body,  Vb,  and  by  the  wake  sheet. 
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Vw  The  velocity  vector,  thus  has  the  form 

V  —  Uoo  +  ^B  +  ^W  —  Uoo  + 

where  (p  is  an  induced  potential  due  to  the  body  and  the  wake. 

In  the  fluid  domain  S,  as  shown  in  Figure  4-1,  the  potential  satisfles  the  Laplace 
equation  ; 

V^(p  =  0.  (4.1) 

Denoting  the  body  surface  as  Sb,  the  kinematic  boundary  condition  on  Sb  is 

V(t)-n  =  -Uoc  n  (4.2) 

Applying  the  Green’s  theorem  with  the  Laplace  equation  (4.1),  the  induced  ve¬ 
locity  can  be  expressed  by  ; 
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(4.3) 


where  q  is  a  variable  point  in  the  integration  and  p  is  a  fixed  point  which  may  be 
located  anywhere  in  the  space.  Sb  is  the  propeller  blade  and  Sw  is  the  wake  surface. 

Equation  (4.3)  shows  that  the  induced  velocity  due  to  the  body  and  the  wake 
can  be  calculated  by  distributing  dipoles  and  sources  on  the  body  and  dipoles  in  the 
wake. 


4.2  Discrete  Formulation 

For  the  numerical  implementation,  the  propeller  surface  and  the  wake  sheet  are  dis¬ 
cretized  into  hyperboloidal  panels  [24].  Constant  strength  sources  and  dipoles  are 
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distributed  on  each  panel  on  the  body  and  bi-quadratic  strength  dipoles  on  each 
panel  in  the  wake.  The  effect  of  the  hub  is  also  included  by  distributing  panels  on 
the  hub  surface.  A  collocation  method  is  applied,  with  the  control  points  being  the 
centroids  of  the  panels.  Using  equations  (3.16)  and  (4.3),  the  discretized  form  of  the 
induced  velocity  becomes  : 
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(4.4) 


where  Nblade  is  the  number  of  blades  and  Npanel  is  the  total  number  of  panels  on  a 
blade,  which  has  N  chordwise  panels  and  M  spanwise  panels  (i.e.  Npanel  =  NxM). 
Nw  is  the  number  of  streamwise  panels  on  each  strip  in  the  wake. 

The  influence  coefficients  aij  and  bi^j  are  deflned  as  the  velocities  induced  at  panel 
7  by  a  unit  strength  dipole  and  source,  respectively,  located  at  panel  j  on  blade  K. 
The  wake  influence  coefficient  IT)^  ^  is  deflned  similarly,  as  the  induced  velocity  at 
panel  i  due  to  the  unit  dipole  at  {m,  1)  in  the  wake  on  the  blade  K.  The  numerical 
computation  of  these  coefficients  are  explained  in  Appendix  C.  The  strength  of  the 
source  (^^)  and  the  dipole  {4>f)  on  the  body  and  the  dipole  {Acf)^,)  in  the  wake 
can  be  found  from  the  typical  panel  method  after  applying  the  kinematic  condition 
on  the  body  and  the  iterative  pressure  Kutta  condition  at  the  trailing  edge  [42],  [24], 
in  which  case  the  geometry  of  the  trailing  wake  sheet  is  assumed  to  be  aligned  with 
the  inflow. 
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4.3  Numerical  Procedure 

In  order  to  include  the  effect  of  the  wake  sheet  roll-up  in  the  flow  adapted  grid,  the 
following  method  is  implemented.  The  method  proceeds  by  computing  the  induced 
velocities  at  the  control  points  in  the  wake,  as  follows. 


Figure  4-2:  The  first  iteration  for  the  FLAG  with  a  wake  sheet  roll-up.  Velocity 
vectors  are  evaluated  at  the  control  points. 

1.  Solve  a  boundary  value  problem  by  using  the  BEM  with  a  geometry  of  the 
trailing  wake  sheet,  which  is  aligned  with  the  inflow.  In  the  case  of  a  wing,  the 
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trailing  wake  will  be  planar  or  cylindrical  (with  very  high  curvature)  and  have 
the  same  direction  as  the  inflow.  In  the  case  of  a  propeller,  the  trailing  wake 
will  be  helicoidal  without  any  contraction  at  the  tip. 


Figure  4-3:  The  second  or  higher  iteration.  Velocity  vectors  are  at  the  control  points. 

2.  Compute  the  induced  velocities  (V^j)  at  the  control  points  of  the  first  row  of 
panels  in  the  wake  from  the  equation  (4.4).  Then  move  the  first  row  of  panels 
to  the  second  row  by  using  a  first  order  Euler  scheme,  which  is  given  as, 

^m,2  ~  ®m,l  iU <x>  T 
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Figure  4-4:  Construct  FLAG  with  wake  sheet  roll-up  from  previous  iteration. 

3.  Adjust  the  downstream  wake  shape  so  that  it  has  the  same  shape  as  the  wake 
at  the  previous  row.  In  other  words,  set  ykj  =  Vij  and  Zkj  —  Zij  for  k  = 
i  -I- 1,  i  +  2,  •  •  • ,  Nw+i,  as  shown  in  Figure  4-2.  Also  check  the  panel  size.  If  the 
panel  size  is  greater  than  the  criteria  Amax,  then  rediscretize  the  panel. 

4.  Repeat  steps  2.  and  3.  until  the  last  row  of  the  panels.  With  this  initial 
geometry,  compute  induced  velocities  at  all  control  points  in  the  wake  and 
move  the  panels  all  together  as  shown  in  Figure  4-3.  Repeat  this  calculation 
until  the  geometry  of  the  wake  sheet  is  converged.  Usually,  two  iterations  have 


been  found  to  be  enough. 


Figure  4-5:  Flow  diagram  for  construction  of  the  FLAG  with  wake  sheet  roll-up  in 
three  dimensions. 

5.  With  the  converged  wake  sheet,  construct  FLAG  and  solve  the  BEM  with  the 
known  wake  sheet  roll-up  surface,  to  get  the  potentials  on  the  body  and  potential 
jumps  in  the  wake,  as  shown  schematically  in  Figure  4-4. 

6.  Repeat  steps  2  to  5  until  the  geometry  of  the  wake  sheet  is  converged. 

The  summary  of  this  procedure  is  shown  in  Figure  4-5.  In  the  calculation  of  the  wake 
sheet  roll-up,  it  usually  takes  two  iterations  with  given  potentials  on  the  body  and 
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the  potential  jumps  in  the  wake.  On  the  other  hand,  for  the  construction  of  FLAG 
with  wake  sheet  roll-up,  two  iterations  are  usually  enough.  As  for  the  CPU,  it  takes 
about  three  times  as  long  as  the  panel  method  without  wake  sheet  roll-up.  In  the 
next  chapter,  this  method  will  be  applied  to  several  wing  and  propeller  geometries 
The  sensitivity  of  the  results  to  the  discretization  parameters  will  be  investigated, 
and  results  of  the  method  will  be  validated  against  existing  experiments. 
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Chapter  5 

Analysis  of  Computational  Results 


In  the  earlier  chapters  of  this  thesis,  the  flow  adapted  grid  and  a  new  wake  sheet  roll¬ 
up  computation  method  were  described  and  incorporated  in  the  panel  method.  In  this 
chapter,  the  new  method  will  be  applied  to  a  number  of  lifting  surface  configurations. 
Two  groups  of  applications  will  be  presented.  The  first  group  is  applications  for  a 
series  of  wings  and  the  second  group  is  applications  for  a  series  of  propellers.  In  order 
to  validate  the  method,  the  results  will  be  compared  to  those  from  other  numerical 
methods  and  to  published  experimental  measurements. 

5.1  Wing 

In  Chapter  3,  a  new  higher  order  panel  method  has  been  developed  in  order  to  predict 
the  geometry  of  the  vortex  sheet  roll-up  in  three  dimensions.  The  method  was  applied 
to  a  lifting  line  with  elliptic  loading.  In  this  section,  in  order  to  examine  the  body 
effect  on  the  geometry  of  the  vortex  sheet  roll-up,  the  method  will  be  applied  to  a 
series  of  wings. 


5.1.1  Rectangular  Wing 

The  method  is  first  applied  to  a  rectangular  wing  of  aspect  ratio  AR=8  at  10  degrees 
angle  of  attack,  as  shown  in  Figure  5-1.  In  this  calculation,  the  panels  in  the  wake 


are  distributed  using  cosine  spacing  in  the  spanwise  and  the  streamwise  directions  to 
emphasize  the  details  near  the  body  and  the  tip.  Following  the  procedure  explained 
in  the  last  chapter,  the  wake  sheet  roll-up  is  obtained  as  shown  in  Figure  5-1. 


INITIAL  GEOMETRY 


AFTER  I-ST  ITERATION 


AFTER  2-ND  ITERATION 


Figure  5-1:  Trailing  wake  sheet  behind  a  rectangular  wing;  AR  =  8,  o;  =  10°, 
{Tlc)max  =  0.01.  40  chordwise  and  30  spanwise  panels  on  the  wing  and  20  streamwise 
panels  in  the  wake.  Constant  thickness  distribution  in  spanwise  direction. 

The  wing  is  discretized  into  40  chordwise  and  20  spanwise  panels.  The  trailing 
wake  sheet  is  discretized  into  20  spanwise  and  20  streamwise  panels.  In  this  particular 
case,  it  took  two  iterations  for  the  wake  shape  to  converge.  The  geometry  of  the  wake 
sheet  at  each  iteration  is  also  shown  in  Figure  5-1.  The  geometries  after  the  first  and 


100 


second  iterations  are  practically  identical  to  each  other.  The  roll-up  pattern  shown 
in  the  bottom  of  Figure  5-2,  has  a  smooth  appearance.  However,  near  the  tip  at  the 
trailing  edge  of  the  wing,  there  is  a  non-smooth  behavior. 


Figure  5-2:  The  geometry  of  the  trailing  wake  sheet  and  its  cross  sections. 

That  may  be  because  the  position  of  the  tip  is  set  at  the  trailing  edge  of  the  tip 
chord.  In  reality,  the  position  of  the  tip  would  be  somewhere  along  the  tip  chord.  The 
present  method  cannot  handle  the  side  edge  separation  yet.  The  geometry  of  the  wake 
sheet  is  compared  to  that  from  Suciu  and  Morino’s  method  [72]  at  two  axial  positions, 
which  are  at  4  and  9  times  the  wing  chord  downstream  from  the  trailing  edge.  Suciu 
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and  Morino’s  method  uses  constant  strength  singularities  and  hyperboloidal  panels 
on  the  wing  and  in  the  wake. 


Z  At  (  X  -  Xj.^  )  1C  -  4.0 
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Figure  5-3:  Wake  cross  section  at  {x  -  xte)  —  4  and  9  x  c  for  a  rectangular  wing  with 
1%  maximum  thickness/chord  ratio,  AR  =  8,  a  10°. 

The  comparison  is  shown  in  Figure  5-3.  Even  though  the  size  of  the  two  roll-up 
regions  is  different,  the  center  of  the  core  and  the  geometry  of  the  rest  of  the  wake 
sheet  shows  good  agreement. 
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5.1.2  An  Elliptic  Wing 

In  Chapter  2,  the  flow  adapted  grid  without  the  wake  sheet  roll-up,  was  applied  to 
several  wings  and  propellers,  and  validated  via  convergence  and  consistency  tests. 


Figure  5-4:  The  geometry  of  the  trailing  wake  sheet  roll-up  behind  an  elliptic  wing 
with  15%.  maximum  thickness/chord  ratio,  a  -  a^eai  =  6°- 

To  examine  the  effect  of  the  vortex  sheet  roll-up  on  the  flow  adapted  grid,  an 
elliptic  wing  is  considered.  The  cross  section  of  the  wing  has  a  NACA662  —  415  shape 
with  an  a=0.8  mean  camber  line.  The  ideal  angle  of  attack  {aideai)  is  -2.5°  and  the 
maximum  thickness/chord  ratio  is  15%.  The  aspect  ratio  is  3.  The  resulting  flow 


adapted  grid  on  the  wing  and  the  geometry  of  the  wake  sheet  roll-up  are  shown  in 
Figure  5-4. 


Figure  5-5:  The  tip  vortex  trajectory  of  an  elliptic  wing  with  15%  maximum  thick¬ 
ness/chord  ratio,  a  —  aideai  —  12.5";  the  thick  line  is  the  tip  vortex  trajectory  from 
the  experiment  given  by  Arndt  (1991). 

The  tip  is  moved  backward  along  the  trailing  edge  and  the  grid  lines  on  the  wing 
seems  to  align  well  with  those  in  the  wake.  In  the  close-up  picture  near  the  tip  of 
the  wing  shown  in  Figure  5-4,  the  tip  vortex  is  shown  to  depart  from  the  tip  of  the 
wing  and  to  roll  up  very  smoothly.  In  this  case,  40  spanwise,  40  chordwise  and  20 
streamwise  panels  are  used  on  the  wing  and  in  the  wake.  It  took  two  iterations  for  the 
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converged  geometry  of  the  wake  sheet  roll-up  and  three  iterations  for  the  flow  adapted 
grid.  Figure  5-5  shows  the  cross  sections  of  the  wake  sheet  roll-up  at  three  different 
locations.  In  the  same  figure,  the  numerical  tip  vortex  trajectory  is  compared  to  that 
from  the  experiment,  which  is  given  by  Arndt  [2]. 


EXPERIMENT 


Figure  5-6:  The  tip  vortex  trajectory  of  an  elliptic  wing  with  15%  maximum  thick¬ 
ness/chord  ratio,  a-(Xideai  =  6.0°  and  12.5°;  the  thick  line  is  the  tip  vortex  trajectory 
from  the  experiment  as  given  by  Arndt(1991). 

The  two  trajectories  agree  very  well  with  each  other.  Arndt  also  showed  that 
in  his  experiment,  the  trajectory  of  the  tip  vortex  did  not  depend  on  the  angle  of 
attack.  In  the  numerical  calculation  of  Krasny  [37]  for  the  same  wing  the  trajectory 


was  found  to  be  dependent  on  the  angle  of  attack,  i.e.  as  the  angle  of  attack  increased, 
the  contraction  of  the  tip  vortex  increased  as  well. 


Figure  5-7:  The  geometry  of  the  trailing  wake  sheet  roll-up  behind  the  swept  elliptic 
wings;  [r/cj^nax  =  0.2,  a  =  5.73^",  aspect  ratio  AR=3.  45°  backward  sweep  (top)  45° 
forward  sweep  (bottom). 

In  the  present  method,  as  the  angle  of  attack  (loading)  decreases,  the  position  of 
the  tip  moves  backward  along  the  wing  trailing  edge.  Thus,  the  tip  vortex  trajectory 
remains  practically  the  same  as  shown  in  Figure  5-6.  In  other  words,  the  loading  does 
not  affect  on  the  tip  vortex  trajectory  as  also  was  observed  in  Arndt’s  experiment. 
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5.1.3  Swept  Elliptic  Wings 

To  examine  the  effect  of  sweep  on  the  geometry  of  the  wake  sheet  roll-up,  the  present 
method  is  applied  on  two  swept  elliptic  wings.  The  first  wing  has  45^*  degree  forward 
sweet  and  the  other  wing  has  45°  backward  sweet. 


Figure  5-8:  The  geometry  of  the  trailing  wake  sheet  roll-up  behind  a  circular  wing 
with  20%  maximum  thickness/chord  ratio,  a  —  O.lrad. 

The  results  for  the  swept  elliptic  wings  are  shown  in  Figure  5-7.  The  aspect  ratio 
of  the  wings  is  3  and  the  angle  of  attack  is  5.73  degrees.  For  both  cases,  the  same 
number  of  panels  are  used  as  that  for  the  non-swept  elliptic  wing.  Since  the  backward 
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swept  wing  has  usually  steeper  slope  in  circulation  distribution  near  the  tip  than  that 
for  the  forward  swept  wing  as  shown  in  Figure  2-13,  the  backward  swept  wing  has  a 
more  advanced  shape  of  the  roll-up  in  the  wake,  as  shown  in  Figure  5-7. 


Figure  5-9:  Circulation  distribution  on  a  circular  wing  planform  hydrofoil;  [r/c]mai  = 
0.2,  a  =  5.73°.  Predicted  by  applying  the  BEM  on  the  flow  adapted  grid  with  roll-up 
;  before  and  after  applying  the  IPK  condition. 


From  the  same  figure,  it  can  be  also  concluded  that  as  the  wing  is  swept  backward, 
the  contraction  angle  of  the  tip  vortex  trajectory  becomes  larger.  This  is  a  different 
trend  from  that  shown  in  Figure  2-15.  It  appears  that  when  the  roll-up  in  the  wake 
is  included  the  final  tip  vortex  contraction  is  smaller  than  that  predicted  when  only 
the  effects  of  thickness  are  included. 
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5.1.4  Circular  Wing 

As  mentioned  in  Chapter  2,  since  a  circular  wing  has  a  very  large  chord  at  the  tip, 
it  is  very  difficult  to  get  a  converged  solution  and  also  to  predict  a  correct  shape  of 
the  wake  sheet  roll-up.  In  particular,  when  an  iterative  Kutta  condition  is  applied  to 
a  thick  circular  wing,  the  circulation  distribution  near  the  tip  becomes  non-physical, 
when  the  conventional  or  the  blade  orthogonal  grid  are  used.  In  this  section,  the 
present  method  is  applied  to  a  circular  wing  with  20%  maximum  thickness/chord 
ratio.  The  angle  of  attack  is  0.1  rad.  The  resulting  shape  of  the  rolled-up  wake  sheet 
is  shown  in  Figure  5-8.  The  position  of  the  tip  moved  more  backward  than  that  in 
the  case  of  the  elliptic  wing.  In  the  same  figure,  the  gridlines  on  the  wing  and  in  the 
wake  are  shown  to  be  smoothly  connected  along  the  trailing  edge  of  the  wing.  The 
tip  vortex  also  leaves  the  tip  and  rolls-up  very  smoothly.  The  circulation  distribution 
is  shown  in  Figure  5-9.  The  difference  between  before  and  after  applying  the  IPK 
condition  is  much  smaller  that  those  from  the  conventional  grid  and  the  flow  adapted 
grid  without  roll-up  as  shown  in  Figures  2-4  and  2-14.  The  flow  adapted  grid  with 
roll-up  has  made  the  necessity  for  the  iterative  pressure  Kutta  condition  less  crucial 
since  the  results  before  and  after  IPK  condition  are  very  close  to  each  other. 

5.2  Propeller 

The  present  method  is  applied  to  propellers  in  the  same  manner  as  in  the  case  of  the 
wings.  The  difference  between  this  computation  and  that  for  the  wing  is  that  each 
trailing  vortex  is  now  assumed  to  travel  in  a  helical  trajectory  rather  than  a  straight 
line  downstream.  Also  for  the  far  field  calculation,  the  sink  disk  is  used  instead  of 
vortex  lines,  which  were  used  in  the  wing  problem.  The  far  field  calculation  starts  at 
X  =  1.5R,  which  is  usually  used  in  the  propeller  application  [16]. 

Consider  a  propeller  subject  to  a  spatially  uniform  flow  'K4(r)  as  shown  in  Figure 
5-10.  The  flow  around  the  propeller  will  be  analyzed  with  respect  to  the  propeller 
fixed  coordinate  system  {x,y,z).  If  the  propeller  rotates  with  angular  velocity  uj, 


z 

Figure  5-10:  Flow  velocity  diagram. 

For  all  the  computations  described  in  this  section,  the  values  of  the  ultimate  wake 
radius  and  the  tip  vortex  contraction  angle  are  computed  by  the  method,  rather  than 
being  given  from  experimental  information. 

Following  the  numerical  procedure  described  earlier,  the  flow  adapted  grid  with 
the  wake  sheet  roll-up,  is  applied  to  three  different  kinds  of  DTMB  propellers.  The 
results  are  compared  to  the  experimental  data. 


no 


5.2.1  Propeller  4990 

The  first  case  is  propeller  4990,  which  is  a  typical  modern  marine  propeller  with  a 
large  hub  and  nonlinear  skew  distribution. 


Figure  5-11:  Trailing  wake  sheet  of  a  propeller  4990. 

In  this  calculation,  only  one  blade  is  considered  and  no  hub  is  included.  The 
propeller  4990  has  a  high  skew  distribution  and  a  big  chord  length  at  the  tip.  The 
calculation  is  performed  at  the  design  advance  coefficient  (J^  =  1.270).  The  number 
of  panel  on  the  blade  is  40  and  20  in  chordwise  and  spanwise  directions,  respectively. 
In  the  wake,  30  streamwise  number  of  panels  are  used.  The  resulting  shape  of  the 
wake  sheet  roll-up  is  shown  in  Figure  5-11.  The  trailing  wake  sheet  leaves  the  trailing 
edge  of  the  blade  and  smoothly  rolls  up  downstream.  In  addition,  the  cross  sections 
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of  the  wake  sheet  are  shown  in  Figure  5-12. 


In  this  figure,  it  is  shown  that  the  core  region  advances  faster  than  the  rest  of 
the  wake  sheet.  It  seems  that  the  axial  flow  in  the  core  becomes  accelerated.  This  is 
explained  in  Appendix  B.1.1. 

5.2.2  Propeller  4119 

The  propeller  4119  is  a  typical  three  bladed  propeller.  The  blade  has  a  big  chord 
length  at  the  tip  and  no  skew  distribution.  The  blade  section  has  an  NACA  66 
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modified  thickness  form  superimposed  on  an  a=0.8  mean  camber  line. 


Figure  5-13;  The  shape  of  the  wake  sheet  of  a  propeller  4119. 


The  calculation  is  performed  at  design  advance  coefficient  equal  to  0.833.  The 
hub  is  not  included  in  this  calculation.  Figure  5-13  shows  a  general  shape  of  wake 
sheet  roll-up  and  also  the  cross  sections.  This  case  also  shows  good  convergence  and 
geometry  of  the  wake  sheet  roll-up.  In  the  cross  section,  unlike  the  propeller  4990,  the 
core  region  is  decelerated.  This  is  also  explained  in  Appendix  B.1.1.  For  this  propeller, 
experimental  results  are  available  in  [27].  Figure  5-14  shows  the  trajectory  of  the  tip 
vortex  along  the  x-direction  for  the  propeller  4990.  The  tip  vortex  trajectory  from 
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the  present  method  gives  an  identical  trajectory  to  that  measured  in  the  experiment. 
In  addition,  the  pitch  of  the  tip  vortex  is  predicted  very  well,  as  shown  in  Figure  5-15. 


r/R 


Figure  5-14:  Radial  location  of  the  tip  vortex  for  a  propeller  4119. 

To  compare  the  circulation  deformations  behind  the  propeller  blades,  two  axial 
locations  are  chosen,  which  are  oX  x/R  =  0.328  and  0.951.  Figure  5-16  shows  these 
two  locations  and  the  circulations  from  the  present  method  and  the  experiment.  The 
present  method  gives  multivalues  of  circulation  distribution  near  the  tip.  That  is 
because  the  circulation  is  calculated  along  the  roll-up  sheet.  In  the  rolled-up  region, 
the  radial  position  of  the  sheet  goes  up  and  down  as  many  times  as  the  number 
of  turns  of  the  roll-up.  On  the  other  hand,  for  the  experiment,  the  circulation  is 
calculated  from  the  distribution  of  mean  tangential  velocity  along  the  cut  in  the 
wake.  This  method  is  described  by  Kerwin  [32]  and  later  Wang  [73]  and  provides 
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a  simple  comparison  between  the  measured  and  the  calculated  circulation.  From 
Kerwin  [32], 

r  _  1  r  Vt 
2'kRV  ~k"r'V 

where  F  is  a  incoming  flow  and  K  is  the  number  of  blades  and  F  is  the  measured 
mean  tangential  velocity. 


Figure  5-15:  Angular  position  of  the  tip  vortex  for  a  propeller  4119. 

The  circulations  from  the  experiment  and  the  calculation  at  two  axial  locations 
agree  well  with  each  other. 
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5.2.3  Propeller  4660 

Finally,  in  order  to  examine  the  hub  effect  on  the  wake  sheet  roll-up,  a  propeller  4660 
is  calculated.  The  propeller  4660  has  five  blades  and  high  skew  distribution. 


Figure  5-16:  Circulation  deformation  in  the  trailing  wake  for  a  propeller  4119;  At 
x/R  =  0.328  and  0.951  . 

The  geometry  of  this  propeller  is  given  in  Table  5.1.  The  calculation  is  performed 
at  J=0.976,  which  is  smaller  than  the  design  advance  coefficient.  In  the  computation, 
the  hub  is  modeled  with  a  fairwater  downstream  and  constant  radius  upstream  as 
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shown  in  Figure  5-17.  In  addition,  the  ultimate  hub  vortex  radius  is  assumed  as 
O.IR.  Figure  5-17  shows  the  predicted  general  shape  of  the  wake  sheet  roll-up.  The 
flow  leaves  the  tip  of  the  blade  and  starts  to  roll  up  smoothly.  Near  the  hub,  the  flow 
follows  the  hub  surface  and  leaves  at  the  ultimate  hub  radius. 


Number  of  Blades 
Hub/Diameter  Ratio 
Section  Meanline 
Section  Thickness  Form 
Design  Advanced  Coefficient 


5 

0.3 

NACA  a=0.8 
NACA66  (Modified) 
1.038 


r/R 

P/D 

^ml  D 

6{degree) 

c/D 

//c 

t/D 

0.30 

1.165 

0.0091 

2.985 

0.178 

0.0000 

0.0420 

0.35 

1.296 

0.0103 

3.481 

0.210 

0.0050 

0.0372 

0.45 

1.480 

0.0103 

4.810 

0.271 

0.0209 

0.0290 

0.55 

1.566 

0.0103 

6.631 

0.327 

0.0267 

0.0226 

0.65 

1.566 

0.0103 

8.978 

0.374 

0.0256 

0.0178 

0.75 

1.498 

0.0103 

11.895 

0.406 

0.0209 

0.0146 

0.85 

1.381 

0.0103 

15.410 

0.409 

0.0151 

0.0122 

0.90 

1.306 

0.0102 

17.403 

0.387 

0.0122 

0.0110 

0.95 

1.222 

0.0103 

19.557 

0.326 

0.0094 

0.0091 

1.00 

1.128 

0.0102 

21.876 

0.000 

0.0000 

0.0000 

Table  5.1:  The  geometry  of  the  propeller  4660. 


From  the  figure,  it  is  clearly  shown  behind  the  propeller  blades  that  as  the  wake 
goes  downstream,  the  tip  vortex  position  is  contracted.  For  this  propeller  experimen¬ 
tal  results  are  available  in  [73].  To  examine  the  change  of  the  circulation  behind  the 
propeller  blades,  three  axial  locations  are  chosen,  which  are  aXx/R  =  0.281,  0.853  and 
1.253.  Figure  5-18  shows  these  three  locations  and  the  circulations  from  the  present 
method  and  the  experiment.  The  agreement  between  the  computed  and  the  mea¬ 
sured  circulation  distributions  behind  the  propeller  blades  is  remarkable.  From  these 
applications,  it  could  be  concluded  that  the  present  method  predicts  the  geometry 
of  the  wake  sheet,  especially  the  tip  vortex  trajectory,  very  smoothly  and  accurately. 
In  addition,  the  computed  circulations  behind  the  propeller  blades  match  well  with 
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Figure  5-17:  The  shape  of  the  wake  sheet  on  propeller  4660. 
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Figure  5-18:  Circulation  deformation  in  the  trailing  wake  for  a  propeller  4660;  At 
x/R  =  0.281,  0.853  and  1.253  . 


Chapter  6 


Conclusions  and 
Recommendations 

In  this  chapter,  the  accomplished  improvements  of  the  propeller  steady  flow  panel 
method  will  be  summarized  and  then  future  research  topics  related  to  the  present 
method  will  be  suggested. 

6.1  Conclusions 

In  this  thesis,  a  new  grid  arrangement,  called  by  “FLow  Adapted  Grid”,  has  been 
developed.  The  flow  adapted  grid  has  been  implemented  in  a  panel  method  and 
applied  to  a  highly  skewed  propeller  and  a  propeller  with  a  large  chord  at  the  tip, 
which  have  had  problems  in  satisfying  the  iterative  pressure  Kutta  condition.  The 
results  from  applying  the  flow  adapted  grid  are  convergent  and  also  consistent  to  the 
results  from  the  lifting  surface  method  [16].  In  particular,  the  3%  discrepancy  in 
circulations  between  the  panel  method  and  the  lifting  surface  method  for  propeller 
4118  with  zero  thickness,  which  was  found  by  Hsin  [24],  is  decreased  to  less  than 
1%.  In  the  case  of  wide  tip  blades  with  no  skew  the  flow  adapted  grid  was  found  to 
improve  the  predicted  pressures  at  the  tip  substantially  and  at  same  time  to  affect 
the  circulation  distribution  appreciabb/.  For  wide  tip  blades  with  high  skew,  the 
proposed  grid  improved  the  predicted  pressures  at  the  tip,  even  though  it  did  not 


affect  the  predicted  circulation  distribution.  In  order  to  validate  the  flow  adapted 
grid,  the  force  coefficients  from  using  the  flow  adapted  grid  are  compared  to  those 
from  existing  experimental  data.  The  resulting  forces  and  tip  vortex  trajectories 
from  using  FLAG  are  found  to  be  closer  to  the  experimental  values  than  those  from 
using  the  conventional  grid.  In  the  flow  adapted  grid  the  position  of  the  tip,  “the 
computational  tip”,  is  determined  by  searching  among  the  streamlines  for  the  one 
which  starts  at  the  largest  blade  radial  location  and  does  not  intersect  the  propeller 
blade.  The  location  of  a  computational  tip  depends  on  the  loading  (angle  of  attack). 

A  robust  and  efficient  three-dimensional  numerical  method  has  been  developed 
for  accurately  predicting  the  geometry  of  the  wake  sheet  roll-up  behind  wings  and 
propellers.  A  hyperboloidal  panel  is  used  in  order  to  improve  the  accuracy  in  the 
highly  roll-up  region,  where  the  curvature  increases.  By  employing  the  hyperboloidal 
panel,  the  gap  between  the  panels  disappears  and  the  surface  of  the  trailing  wake 
sheet  can  be  modeled  more  accurately.  In  addition,  in  order  to  avoid  discontinuities 
in  the  dipole  strength  in  the  wake,  bi-quadratic  strength  dipoles  are  used.  The  panel 
method  can  predict  the  shape  of  the  wake  sheet  smoothly.  For  the  more  advanced 
shape  of  the  wake,  a  rediscretization  scheme  is  applied  at  each  iteration.  In  the  first 
iteration,  the  vortex  sheet  moves  from  row  to  row  in  a  stripwise  sense.  From  the 
second  and  higher  iteration,  the  vortex  sheet  moves  altogether.  The  results  of  the 
high-order  panel  method  are  convergent  and  predict  a  smooth  and  reliable  geometry 
of  the  wake  sheet  roll-up.  The  results  of  the  method  are  also  validated  to  those  from 
other  numerical  methods. 

Finally,  the  effect  of  the  three-dimensional  vortex  sheet  roll-up  is  included  in 
the  flow  adapted  grid.  This  method  uses  an  iterative  scheme  to  satisfy  the  boundary 
condition  in  the  wake,  where  the  pressure  jump  is  forced  to  be  equal  to  zero.  For  each 
iteration,  the  geometry  of  the  wake  sheet  roll-up  is  obtained  with  a  given  circulation 
distribution,  which  is  then  updated  by  the  panel  method  applied  on  the  blade  with  the 
rolled-up  wake  sheet.  The  results  of  this  scheme  show  good  convergence  with  number 
of  panels.  In  addition  it  has  been  found  that  when  the  FLAG  and  the  completely 
rolled-up  wake  geometry  are  combined  the  results  before  and  after  IPK  condition 
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become  very  close  to  each  other.  The  method  is  applied  to  many  wing  and  propeller 
geometries  and  the  results  are  shown  to  agree  with  the  existing  experimental  data. 


Figure  6-1:  Incorrect  location  of  the  tip. 

6.2  Recommendations  for  Future  Research 

The  following  improvements  on  both  the  physical  modeling  and  the  numerical  scheme 
of  the  present  method  can  be  made. 
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®  When  a  tip  is  assumed  to  be  at  the  highest  radial  position,  the  wake  sheet  will 
not  be  contracted  but  instead,  it  will  diverge  as  shown  in  Figure  6-1.  From  this 
examination,  the  position  of  the  tip  vortex  detachment  point,  “the  computa¬ 
tional  tip”,  may  be  determined  numerically  as  shown  in  Figure  6-2. 


Figure  6-2:  The  tip  vortex  detachment  point. 

Assuming  a  position  of  the  tip,  the  total  velocity  {V total)  is  computed  at  a 
control  point  of  the  panel  in  the  wake,  which  is  just  next  to  the  assumed  tip. 
If  the  radial  component  of  the  total  velocity  points  outwards  then  move  the  tip 
downstream  along  the  trailing  edge  of  the  blade  until  the  rolled-up  wake  at  the 
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tip  leaves  the  blade  smoothly.  This  scheme  will  work  well  but  the  panel  method 
will  have  to  be  solved  for  each  assumed  tip.  This  will  increase  the  CPU-times 
appreciably. 


Figure  6-3:  Initial  paneling  for  the  leading  edge  separation  for  an  elliptic  wing. 

•  The  leading  edge  separation  may  be  included  by  moving  the  tip  forward  along 
the  leading  edge  of  the  propeller  blade.  The  initial  panel  arrangement  is  shown 
in  Figure  6-2.  Applying  IPK  condition  at  the  trailing  edge  and  the  part  of 
leading  edge  downstream  of  the  detachment  point,  the  circulation  distribution 
can  be  computed.  With  this  circulation,  the  geometry  of  the  trailing  wake  sheet 
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and  the  wake  sheet  separated  from  the  leading  edge  may  be  computed. 

Other  improvements  which  are  recommended  for  further  research  are  : 

®  The  modeling  of  a  viscous  core  and  the  calculation  of  the  boundary  layer  thick¬ 
ness  at  the  tip. 

•  The  prediction  of  the  leading  edge  separation  start  point. 

®  The  modeling  of  tip  vortex  cavitation. 
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Appendix  A 


A  Numerical  Kutta  Condition 

The  circulation  around  the  body,  which  is  equal  to  the  potential  jump  in  the  wake,  is 
determined  by  applying  the  Kutta  condition  at  the  trailing  edge.  The  Kutta  condition 
requires  that  the  velocity  at  the  trailing  edge  be  finite.  In  the  numerical  calculation, 
the  first  approximate  Kutta  condition,  which  is  called  Morino’s  condition  [59],  requires 
that  the  strength  of  the  dipole  in  the  wake  be  equal  to  the  difference  in  the  value  of  the 
dipole  strengths  of  the  two  panels  at  the  trailing  edge.  However,  this  condition  does 
not  include  the  free  stream  (cross  flow),  which  is  in  the  direction  of  a  line  connecting 
the  control  points  of  the  two  panels  at  the  trailing  edge.  An  iterative  pressure  Kutta 
condition  (IPK  condition)  is  thus  introduced  with  the  free  stream  correction  [42],  [24]. 
The  strength  of  the  dipole  in  the  wake  is  adjusted  until  the  pressure  on  the  upper 
and  lower  panels  become  equal  at  the  the  trailing  edge.  However,  some  problems 
were  found  when  this  condition  was  applied  to  analyze  a  high  skewed  propeller  or  a 
propeller  with  a  large  tip  chord  [24].  In  this  appendix,  these  problems  will  be  discussed 
and  improvements  to  the  numerical  Kutta  condition  will  be  presented.  These  issues 
were  first  addressed  by  Kinnas  et  al.  [36]  but  will  be  repeated  in  this  Appendix  for 
the  sake  of  completeness. 

The  upper  and  lower  panels  at  the  trailing  edge  are  considered.  The  total  velocity 
vectors  on  both  sides  of  the  trailing  edge  are  given  as  and  V~  at  the  suction 
and  pressure  sides  of  the  trailing  edge,  respectively,  as  shown  in  Figure  A-1.  The 
direction  of  the  vorticity  vector  7  in  the  wake  is  also  shown  in  the  same  figure.  The 
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two  velocity  vectors  have  equal  magnitudes  as  a  result  of  the  IPK  condition^: 


V 


+ 


V- 


(A.l) 


On  the  other  hand,  the  direction  of  the  mean  velocity  vector,  Vm  —  [V^  +  V  ]/2, 
is  very  different  from  that  of  the  wake  vorticity  vector  7. 


Figure  A-1:  Velocity  vectors  on  the  suction  and  pressure  sides  at  the  trailing  edge  of 
a  circular  planform  wing;  [r/cj^oi  =  0.2,  a  =  O.lrad.  After  an  IPK  condition. 

This  will  result  in  a  pressure  jump  in  the  wake  given  by: 


Apw  =  p\V,n  X  7l 


(A.2) 


^The  IPK  condition  was  found  to  affect  the  magnitude  of  the  trailing  edge  velocities  more  than 
their  directions 
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where  p  is  the  flow  density. 

In  other  words,  even  though  the  IPK  condition  has  ensured  the  equality  of  pres¬ 
sures  at  the  trailing  edge  on  the  blade,  it  has  no  way  to  force  the  zero  pressure  jump 
condition  in  the  wake.  Instead,  we  have  to  align  the  vorticity  vector,  i.e.  wake  geom¬ 
etry,  with  the  mean  velocity  vector  at  the  trailing  edge.  It  would  also  seem  natural 
for  the  grid  on  the  blade  to  be  aligned  with  the  mean  velocity  vector  at  the  trailing 
edge.  This  can  be  explained  as  follows. 


Figure  A-2;  Schematic  of  the  paneling  on  a  3-D  hydrofoil  and  its  wake  in  the  vicinity 
of  the  trailing  edge. 

1 
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The  total  potentials  and  $  at  the  suction  and  pressure  sides  at  the  trailing 
edge  ,  respectively,  may  be  expressed  as  follows: 

$+  =  + 

$-  =  ^^  +  V-As  (A.3) 

where  and  are  the  total  potentials  at  the  control  points  of  the  trailing  edge 
panels  at  the  suction  and  pressure  sides,  respectively;  As  is  the  distance  of  the  control 
points  from  the  trailing  edge  measured  along  the  “chordwise”  grid  direction  on  the 
planform,s,  as  shown  in  Figure  A-2;  and  V~  are  the  projections  of  the  total 
trailing  edge  velocities  V'*'  and  V~,  respectively,  along  s. 

The  Morino  Kutta  condition  [59]  at  the  trailing  edge  is: 

A<f>  =  $+  -  (A.4) 

On  the  other  hand,  the  numerical  implementation  of  equation  (A.4)  requires: 

A$m  =  $1  —  "Fyv  (A.5) 

Thus,  the  discretization  error,  E,  in  implementing  the  Morino  condition  (i.e.  before 
the  IPK  condition),  may  be  expressed,  by  making  use  of  equations  (A.3),  (A.4)  and 
(A.5),  as  follows: 

^  =  A$  -  A$m  =  As  [F+  -  1/-]  (A.6) 

According  to  equation  (A.6),  in  order  to  minimize  the  difference  between  the  circu¬ 
lations  before  and  after  applying  the  IPK  condition,  we  should  have: 

=  W  (A-7) 

In  light  of  equation  (A.l),  equation  (A. 7)  is  equivalent  to  requiring  that  the  mean 
velocity  vector  Vm  is  aligned  with  the  grid  direction  s  on  the  planform.  This  explains 
the  larger  difference  between  the  circulation  distributions  before  and  after  the  IPK 
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condition,  in  the  case  of  the  blade  orthogonal  grid  than  in  the  case  of  the  conventional 
grid,  stated  earlier.  In  the  case  of  the  blade  orthogonal  grid  the  angle  between  Vm 
and  the  $  direction  is  larger  than  in  the  case  of  the  conventional  grid. 
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Appendix  B 


Modeling  of  The  Roll-Up  Region 

This  appendix  is  a  brief  summary  of  Lee  [43]  and  Cummings  [8]  work,  which  explains 
the  rationale  behind  modeling  the  roll-up  region  as  a  core,  a  subcore  and  a  sheath. 
The  nature  of  each  region  will  be  examined  and  the  pressure  calculation  in  the  core 
will  be  expressed. 

B.l  Definition 

The  clue  toward  a  description  of  the  roll-up  region  is  provided  by  experiments  [10].  In 
this  section,  these  three  regions,  referred  to  as  the  core,  the  subcore  and  the  sheath, 
are  defined  as  follows  : 

•  Core  ;  small  changes  in  the  tangential  velocity. 

•  Subcore  ;  large  changes  in  the  tangential  velocity  close  to  the  vortex  axis. 

•  Sheath  ;  moderate  changes  in  the  tangential  velocity  surrounding  the  core. 

These  three  regions  are  the  outcome  of  the  vortex  sheet  roll-up.  The  following 
summarizes  each  region  : 
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B.1.1  The  Core 


As  shown  in  Figure  B-1,  an  orthogonal  curvilinear  coordinate  system  can  be  estab¬ 
lished  in  which  the  core  axis  s  coincides  with  one  of  principal  axes  of  this  curvilinear 
system.  The  other  coordinates  (r,  0)  are  defined  in  a  plane  normal  to  the  local  tangent 
vector  to  s. 


r 


Figure  B-1:  The  coordinate  system  {s,r,9)  along  the  core  axis(s) 


Then,  the  non-dimensionalized  Navier-Stokes  equation  for  steady  incompressive 
laminar  flow  in  this  coordinate  system  are  ; 


Vr 


dVr 

dr 


-f 


4- 


Ve  dVr  Vg  dVr  _  dp 

r  do  r  ds  dr 

1  d  il  d  .  \  1  d'^Vr  d^Vr 

Re  dr  yr5r  r'^  dO'^  ds"^ 


2  dvg 

r^  do 


(B.l) 
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dve  vedve  VrVe  dve 
dr  r  09  r  ^  ds  r  06 

1  \  0  f  1  0  ,  A  1  d'^ve  Q'^ve  2  dvr 
Re  dr  \r  dr  ^  y  09"^  ds^  09 


(B-2) 


Vr 


dvs 

dr 


vedvs  , 

dvs  _ 

-b 

r,a  + 

r  09 

ds 

+ 

1 

dvg 

Re 

r  dr  \ 

^  dr 

dp 

ds 

1  d'^Vs 

+  J.2  QQ2 


d^Vs 

ds^ 


(B.3) 


The  core  is  defined  as  the  region  with  a  relatively  flat  tangential  velocity  profile. 
From  this  definition,  gradients  in  the  tangential  velocity  are  small.  Also  a  Reynolds 
number  can  be  assumed  to  be  large.  Therefore,  at  equation  (B.2),  the  nondimen- 
sional  laminar  viscous  term  (r becomes  much  smaller  that  the  convec¬ 
tive  terms.  In  addition,  measurements  from  Schmucker  and  Gersten  [70]  indicate 
that  in  the  core  region  the  mean  square  quantities  of  the  turbulent  fluctuations  are 
estimated  as  about  3%.  It  can  be  concluded  that  both  laminar  and  turbulent  shear 
stresses  are  small  in  the  core,  so  that  the  core  can  be  regarded  as  inviscid. 

Inside  the  core,  the  axial  velocity  is  usually  accelerated  [38],  [44].  A  simple  expla¬ 
nation  of  this  is  follows.  In  Figure  B-2,  at  two  axial  stations,  A  and  B,  the  pressures 
on  the  same  streamline,  have  the  same  value  if  the  flow  field  in  the  core  is  assumed 
conical. 

P2A  =  P2B 

,where  the  first  subscript  refers  to  the  position  on  the  same  streamline,  and  the  second 
indicates  the  axial  station.  The  radial  momentum  equation  is 


dp 
dr 
1  dp 
pdr 


=  p{v  ■  Vv)r 


Vr~ 


1  dVr 

r  09 


-  —  +  v, 

r 


(B.4) 
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Assuming  a  slender  vortex  core,  the  radial  momentum  equation  (B.4)  reduces  to 


dp 

dr 


>  0 


Thus 


P2A  >  Pi  A 


is  obtained.  This  implies  that  there  exists  a  favorable  axial  pressure  gradient,  which 
accelerates  the  flow  in  the  core. 


Figure  B-2:  Axial  flow  in  the  Core 

This  happens  only  if  the  conical  approximation  and  the  slender  core  assumption 
are  valid.  For  the  propeller  4990,  the  flow  near  the  tip  can  be  assumed  as  a  conical 
flow  because  the  propeller  blade  has  high  skew  distribution  and  large  tip  chord.  In 
addition  the  trailing  wake  sheet  is  not  much  rolled-up  as  shown  in  Min’s  experiment 
[54],  in  which  the  contraction  is  decreased  as  the  skew  becomes  larger.  That  may  be 
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a  reason  why  the  flow  in  the  core  region  is  accelerating  in  the  numerical  results  as 
shown  in  Figure  5-11.  On  the  other  hand,  propeller  4119  does  not  have  skew.  The 
flow  near  the  tip  is  similar  to  that  around  the  circular  wing.  In  this  case,  the  conical 
flow  approximation  is  not  valid  anymore.  Also  since  the  trailing  wake  sheet  has  a 
bigger  core  size  than  that  for  the  propeller  4990,  the  slender  core  assumption  may 
not  be  valid.  That  may  be  one  reason  that  in  this  case  the  flow  is  decelerated  in  the 
core,  as  shown  in  Figure  5-12. 

B.1.2  The  Subcore 

The  subcore  is  a  small  region,  which  is  deep  inside  the  core  and  close  to  the  vortex 
axis.  In  this  region,  the  tangential  velocity  changes  rapidly.  Inside  the  subcore, 
viscous  effects  can  be  no  longer  neglected  because  laminar  viscous  terms  become 
large  as  r  decreases.  In  this  thesis,  this  region  is  ignored. 


BOUNDARY  LAYER  THICKNESS 


Figure  B-3:  Schematic  diagram  showing  a  roll-up  shear  layer  and  the  cut-off  location 


135 


B.1.3  The  Sheath 


This  is  the  region  shielding  the  core  from  the  external  flow.  Unlike  the  core  the 
sheath  is  a  region  in  which  both  gradients  in  tangential  velocity  and  viscous  effects 
are  substantial.  The  sheath  itself  is  just  a  remnant  of  the  wake  sheet  in  the  roll-up 
process.  However,  physically,  this  is  the  only  medium  from  which  irrotational  flow, 
entering  the  rotational  but  inviscid  core,  can  receive  vorticity,  thereby  preserving 
Kelvin’s  theorem.  In  other  words,  if  the  sheath  does  not  exist,  the  solutions  inside 
the  core  and  outside  the  core,  cannot  be  matched. 


B.2  The  Size  of  the  Core 

This  section  will  explain  how  the  size  of  a  core  is  determined.  The  center  of  the  core 
coincides  with  the  position  of  the  trajectory  of  the  tip  vortex,  which  is  determined  by 
finding  the  centroid  of  the  points  near  the  core  in  the  numerical  calculation.  The  core 
boundary  in  a  plane  normal  to  the  tip  vortex  trajectory  is  approximated  as  a  circle. 
When  the  boundary  layer  thickness  is  added  to  the  wake  sheet  at  the  tip,  there  will 
be  a  junction  location  at  which  neighboring  arcs  begin  to  touch  ,Co  in  the  Figure  B-3. 
The  core  radius  is  defined  as  the  distance  between  the  core  center  and  the  junction 
point.  Since  the  viscous  calculation  around  the  vortex  core  is  beyond  the  limit  of  the 
present  thesis,  this  will  not  be  implemented  in  this  thesis. 


B.3  The  Pressure  in  the  Subcore 


The  accurate  calculation  of  the  pressure  coefficient  in  the  core  is  essential  in  the 
prediction  of  the  tip  vortex  cavitation  inception  because  the  minimum  pressure  is 
inside  the  core.  In  this  section,  the  calculation  of  the  pressure  inside  the  core  will  be 
explained. 

Assume  the  motion  is  axisymmetric.  The  equation  of  motion  follows. 


u 


du 

dx 


du 


1  dp 
pdx 


-f-  u'V'^u 
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dv  dv  nP' 

VT 

dw  dw  w 

- ^  - ^ 

ax  or  r 


=  p 


1  dp 


V^-y- 


ti; 


(B.5) 


Assume  axial  gradients  are  small  compared  to  radial  gradients.  Thus, 

A  A 

dx  dr 

u  «  V 

Then  equations  (B.5)  become, 

(B.6) 
(B.7) 
(B.8) 

The  integral  form  of  equation  (B.7)  is  ; 

TOO  y;2 
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du  du  1  dp 
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In  the  region  far  downstream.  Bachelor  finds  a  solution  for  this  equation  of  form: 


w  — 


r 

2'Kr 


l-l 


If  the  subcore  radius  is  defined  as  the  radius  where  the  tangential  velocity  w  is  max¬ 
imum,  then  the  subcore  radius  a  has  a  following  form. 


a 


1.258 


Aux 


The  change  of  subcore  radius  with  the  distance  downstream  is  ; 

da  a 
dx  2x 
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From  the  above  results,  it  can  be  concluded  that  the  subcore  radius  does  not 
change  fast  enough  for  viscosity  to  have  any  significient  effect  on  results  in  the  region 
close  to  the  propeller  blade.  Therefore,  the  subcore  can  be  modeled  by  a  Rankine 
vortex  with  a  group  of  trailing  vorticies  outside  it,  at  each  downstream  section.  The 
pressure  in  the  subcore  can  be  calculated  assuming  axial  symmetry  and  a  constant 
radius  sub  core. 

For  r  <  a 


p{r)  =  Pa,- P 
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Therefore,  the  pressure  coefficient  Cp  is  ; 
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Appendix  C 


The  Calculation  of  The  Induced 
Velocity 


In  the  calculation  of  the  three-dimensional  trailing  wake  sheet  roll-up,  a  new  higher 
order  panel  method  is  applied.  This  method  is  developed  in  order  to  handle  the 
arbitrary  strength  of  singularities  over  a  curved  surface.  A  special  algorithms  has  been 
developed  for  the  calculation  of  induced  potential  by  Maniar  [45].  In  this  appendix, 
the  induced  velocity  due  to  a  distribution  of  normal  dipoles  on  a  curved  panel  is 
evaluated  for  the  case  where  the  density  of  the  singularities  is  of  arbitrary  polynomial 
form.  The  same  algorithm  as  Maniar  [45]  will  be  applied  for  the  calculation. 

As  shown  in  Figure  C-1,  an  orthogonal  coordinate  system  {x,y,z)  is  taken  as  a 
global  coordinate  system.  Then  a  panel  is  transformed  to  local  coordinate  system 
(^)  Vi  C)  with  C  =  0.  Let  S  be  the  curved  panel  surface  and  S  be  its  transformed  panel 
surface.  Q  is  a  point  on  S  and  P  is  a  field  point.  Then  the  induced  velocity  at  P  due 
to  a  distribution  of  dipoles  of  order  M  on  S  is  given  as 


(C.l) 


Figure  C-1:  The  panel  geometries  in  the  global  coordinate  system  and  the  transformed 
coordinate  system. 

Consider  the  evaluation  of  gradient  of  a  weighted  dipole  integral  in  (C.l), 
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The  equation  (C.2)  can  be  rewritten  in  the  (f,  t?,  C)  plane  as, 


(C.2) 


r  j  ,■  r  nR^  -  3R{n  ■  R) 


P  H„7  r  o-^  •  (^  •  (^C  ^  ^r,)) 

L  4 


d^dr]  (C.3) 


This  integral  will  be  implemented  in  three  categories,  which  are  near  field,  far 
field  and  self  induced  coefficients.  Each  region  is  decided  by  the  ratio  between  the 
maximum  diagonal  distance  of  the  panel  and  the  distance  between  the  field  point 
and  the  collocation  point  of  the  panel.  The  followings  summarize  the  numerical 
formulation  for  each  region; 


•  Far  Field 

This  is  the  case  when  the  field  point  is  far  away  from  the  panel.  The  ratio  of 
the  distance  between  the  field  point  and  the  collocation  point  to  the  maximum 
diagonal  distance  is  greater  than  5.  Usually  all  panels  except  themselves  and 
neighboring  panels  are  in  this  category.  The  integrals  in  equation  (C.3)  do 
not  have  a  closed  form  solution.  Therefore,  the  integrals  have  to  be  solved 
numerically.  Gaussian  Quadrature  is  used  in  this  section. 

The  discrete  formulations  are  follows; 


C  cii  (xg  X  _  qR'  [R  '  ^ 
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6-6\  fV2-vi 


where  ^ 


6  ~  6  ,  6  +  6 

V2-V1  /  ,  ^1+^2 
— 2 — ^  - 2 — 

x„  —  X 


is  the  order  of  quadrature, 
are  coefficients  of  the  quadrature. 


In  equation  (C.4),  the  directed  surface  element  x^  x  Xr,  can  be  computed  by 
follows. 
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2N-1  min(n,N) 
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^=0 

From  equations  (C.4)  and  (C.5),  the  induced  velocity  can  be  computed. 

•  Near  Field 

In  near  field,  the  induced  velocity  is  evaluated  by  an  adaptive  subdivision 
scheme,  which  will  be  explained  below.  When  the  field  point  is  in  the  near 
field  region,  the  panel  is  subdivided  until  the  ratio  between  the  distance  from 
the  field  point  to  the  collocation  point  and  the  maximum  diagonal  distance  of 
the  subdivided  panel  is  out  of  the  near  field  region.  Then,  we  can  use  the  same 
scheme  for  that  for  the  far  field.  Let  the  number  of  subdivided  panels  along  the 
^  and  7]-direction  be  N^,Njj  respectively.  Then 

_ ^  _ II  rVn^m-i-l  f^n+ly7n  . _ 

=  £X  /  /  (C.6) 

71=1  m=l  •^?n,m  -ti 

When  the  panel  is  subdivided,  and  Nr/  are  determined  from  the  panel  aspect 
ratio.  The  Equation  (C.6)  can  be  computed  by  using  the  Gaussian  quadrature 
as  explained  in  case  for  the  far  field  calculation. 

•  Self  Induced  Coefficient 

This  is  the  case  when  the  field  point  is  on  the  surface  over  which  integration  is 
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carried  out.  For  the  induced  potential,  the  singularities  are  analytically  removed 
and  the  desingularized  integrands  are  expanded  in  a  series  with  algebraic  type 
terms  which  are  integrable  [45]. 


Z 


Figure  C-2:  An  example  for  calculating  self  induced  velocity. 

In  the  case  of  self  induced  velocity,  to  this  moment,  there  has  been  found  no  way 
to  remove  the  singularities  analytically.  Instead,  we  divide  a  panel  into  smaller 
panels  with  the  center  panel  having  the  same  collocation  point  as  the  original 
panel.  Then,  find  the  self  induced  velocity  by  assuming  that  the  strength  of  the 
dipole  is  constant  on  the  central  panel  and  its  geometry  flat.  In  the  other  panels 
the  induced  velocities  are  computed  by  using  the  same  scheme  as  that  for  the 
near  and  the  far  field.  In  order  to  find  an  adequate  number  of  subdivisions,  an 
example  is  considered  for  a  non-planar  panel  with  bi-linear  dipole  distribution 
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as  shown  in  Figure  C-2.  As  shown  in  Table  C.l,  the  self  induced  velocity  has 
converged  in  3  subdivisions.  In  the  calculation  of  the  wake  sheet  roll-up,  three 
subdivision  is  used  for  the  shorter  side  and  for  the  longer  side,  the  number  of 
subdivisions  is  determined  by  the  aspect  ratio  of  a  panel. 


Number  of  Subdivisions 

u 

V 

w 

1 

0.841592 

0.841592 

-1.683183 

3 

0.819042 

0.810093 

-1.642820 

5 

0.819042 

0.810093 

-1.642820 

Table  C.l:  Convergence  of  the  self  induced  velocity. 
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